Prepare the workspace: clean, set options, load packages, load and final formatting on inputs.
rm(list = ls()) # clear workspace
knitr::opts_chunk$set(tidy = TRUE, collapse = TRUE) # set global knitr options.
# load required packages
require(ggplot2); # used for plotting
require(pwr); # used for power analysis
require(lme4); # used for mixed effects models
require(lmerTest); # used to interpret mixed effects models
require(sjPlot); # used to output tables
require(weights); # used for weighted t.test
require(plyr); # used for dataframe manipulation
options(scipen = 999) # turn off scientific notation
getwd()
## [1] "/Users/boterolab1/Box Sync/CB_VF_Shared/Wet_Lab/Projects/MP_Clean/Methods_Manuscript_V6/File_S3/3_Analysis"
dir_in_pocpro <- "2_CountstoAnalysis/001_Prepare_Datasets" # input directory for Proof of concept processed dataframes
setwd("..")
setwd(dir_in_pocpro) # set to above directory path
load(file = "poc1.Rdata") # load the 1bc control processed dataframe
load(file = "poc2.Rdata") # load the 2 bc control processed dataframe
load(file = "poc92.Rdata") # load the 92 bc proof of concept fitness assays processed dataframe
dir_in_pocraw <- "001_Raw_Data_Metadata" # input directory for Proof of concept raw dataframes
setwd("..")
setwd(dir_in_pocraw) # set to above directory path
load(file = "001_counts.expected.Rdata")
poc.ce <- (counts.expected + 1)[c(1:92), ]
rm(counts.expected) # expected counts
load(file = "001_counts.unexpected.Rdata")
poc.cc <- (counts.unexpected + 1)[c(1:92), ]
rm(counts.unexpected) # cross contamination counts.
dir_in_faevoraw <- "002_Raw_Data_Metadata" # input directory for FA and Evo raw dataframes
setwd("..")
setwd(dir_in_faevoraw) # set to above directory path
load(file = "002_CC_counts.expected.Rdata")
faevo.ce <- (ects + 0)
rm(ects) # expected counts {combined}
load(file = "002_CC_counts.unexpected.Rdata")
faevo.cc <- (ccts + 0)
rm(ccts) # cross contamination counts.{combined}
load(file = "002_counts.expected.Rdata")
faevo.1.ce <- (counts.expected + 1)[c(1:78), ]
rm(counts.expected) # expected counts {first run samples}
load(file = "002_counts.unexpected.Rdata")
faevo.1.cc <- (counts.unexpected + 1)[c(1:78), ]
rm(counts.unexpected) # cross contamination counts. {first run samples}
load(file = "002_Reruns_counts.expected.Rdata")
faevo.2.ce <- (counts.expected + 1)[c(1:78), ]
rm(counts.expected) # expected counts {rerun samples}
load(file = "002_Reruns_counts.unexpected.Rdata")
faevo.2.cc <- (counts.unexpected + 1)[c(1:78), ]
rm(counts.unexpected) # cross contamination counts. {rerun samples}
dir_in_faevopro <- "002_Prepare_Datasets" # input directory for FA and Evo processed dataframes
setwd("..")
setwd(dir_in_faevopro) # set to above directory path
load(file = "myfa.Rdata") # fitness assay
load(file = "myevo.Rdata") # evolution
rm(dir_in_faevopro, dir_in_faevoraw, dir_in_pocpro, dir_in_pocraw)
Create helper functions for linear model residual plotting and data transformation.
# expand_residuals creates a modified histogram of residuals. Each entry in
# the model output is plotted a number of times equal to its reads.
expand_residuals <- function(model, reads) {
myfr <- NA
for (i in 1:length(residuals(model))) {
times <- (((reads[i] - min(reads, na.rm = T))/(max(reads) - min(reads))) *
1000) + 1
myfr <- c(myfr, rep(residuals(model)[i], times = times))
}
return(myfr)
}
# log transformation for modeling.
log_transform <- function(myvar) {
log(((myvar - min(myvar))/(max(myvar) - min(myvar))) + 1)
return(myvar)
}
Set color palette options for all plots
setpalette <- RColorBrewer::brewer.pal(6, "Dark2") # current preferred palette
# setpalette <-
# c('#993404','#d95f0e','#fe9929','#41b6c4','#2c7fb8','#253494') # option 2
# setpalette <- c('#8c510a','#d8b365','#f6e8c3',
# '#c7eae5','#5ab4ac','#01665e') # option 3
Analysis: Assess POC fitness assay barcode fitness equivalence.
my <- poc92[order(poc92$bcid), ]
# model & diagnostic plots
# ----------------------------------------------------
my$rw_corrected <- my$rw - weighted.mean(my$rw, my$r)
mymod <- lm(my$rw_corrected ~ my$bcid + 0, weights = my$r)
summary(mymod)
##
## Call:
## lm(formula = my$rw_corrected ~ my$bcid + 0, weights = my$r)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -4.7350 -0.5858 -0.0109 0.5304 5.3405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## my$bcidd1A10 0.00132555 0.00371437 0.357 0.721280
## my$bcidd1A11 -0.00762370 0.00466055 -1.636 0.102267
## my$bcidd1A2 0.00811450 0.00396019 2.049 0.040778 *
## my$bcidd1A3 0.00494242 0.00439553 1.124 0.261165
## my$bcidd1A4 0.00795573 0.00448555 1.774 0.076495 .
## my$bcidd1A5 0.00103320 0.00520070 0.199 0.842573
## my$bcidd1A6 -0.00085773 0.00357809 -0.240 0.810610
## my$bcidd1A7 0.00344369 0.00512097 0.672 0.501476
## my$bcidd1A8 0.00583557 0.00414724 1.407 0.159778
## my$bcidd1A9 0.00724382 0.00374750 1.933 0.053584 .
## my$bcidd1B1 -0.00527477 0.00502129 -1.050 0.293807
## my$bcidd1B10 -0.00821202 0.00615734 -1.334 0.182674
## my$bcidd1B11 -0.01330702 0.00511428 -2.602 0.009437 **
## my$bcidd1B12 -0.00741678 0.00553232 -1.341 0.180413
## my$bcidd1B2 0.00289719 0.00431315 0.672 0.501957
## my$bcidd1B3 0.00354456 0.00434513 0.816 0.414878
## my$bcidd1B4 0.00951091 0.00421461 2.257 0.024293 *
## my$bcidd1B5 0.00961576 0.00421606 2.281 0.022819 *
## my$bcidd1B6 0.00134667 0.00501745 0.268 0.788462
## my$bcidd1B7 -0.00398474 0.00435334 -0.915 0.360287
## my$bcidd1B8 0.00386754 0.00435584 0.888 0.374856
## my$bcidd1B9 -0.00197406 0.00517408 -0.382 0.702909
## my$bcidd1C1 -0.00130215 0.00523917 -0.249 0.803778
## my$bcidd1C10 -0.00200878 0.00612990 -0.328 0.743221
## my$bcidd1C11 0.00332553 0.00513862 0.647 0.517707
## my$bcidd1C12 -0.00495638 0.00525931 -0.942 0.346266
## my$bcidd1C2 0.00310014 0.00440966 0.703 0.482235
## my$bcidd1C3 0.00170870 0.00544151 0.314 0.753591
## my$bcidd1C4 0.00064326 0.00445091 0.145 0.885124
## my$bcidd1C5 -0.00293377 0.00470909 -0.623 0.533457
## my$bcidd1C6 -0.00108375 0.00401081 -0.270 0.787069
## my$bcidd1C7 -0.00224296 0.00857907 -0.261 0.793814
## my$bcidd1C8 -0.00815331 0.00535937 -1.521 0.128566
## my$bcidd1C9 0.01983075 0.00598788 3.312 0.000968 ***
## my$bcidd1D1 0.00371456 0.00415556 0.894 0.371651
## my$bcidd1D10 0.01036458 0.00405675 2.555 0.010802 *
## my$bcidd1D11 -0.00004981 0.00380942 -0.013 0.989570
## my$bcidd1D12 0.00301900 0.00386910 0.780 0.435449
## my$bcidd1D2 -0.01098378 0.02236414 -0.491 0.623464
## my$bcidd1D3 -0.00440055 0.00562205 -0.783 0.434011
## my$bcidd1D4 -0.00693858 0.00737022 -0.941 0.346760
## my$bcidd1D5 -0.00326940 0.00563709 -0.580 0.562087
## my$bcidd1D6 -0.00844038 0.00554061 -1.523 0.128053
## my$bcidd1D7 -0.00168313 0.00462568 -0.364 0.716052
## my$bcidd1D8 0.00180750 0.00397535 0.455 0.649462
## my$bcidd1D9 0.00014080 0.00545635 0.026 0.979419
## my$bcidd1E1 -0.01270595 0.00523534 -2.427 0.015441 *
## my$bcidd1E10 -0.00147548 0.00433537 -0.340 0.733691
## my$bcidd1E11 -0.00708010 0.00406068 -1.744 0.081609 .
## my$bcidd1E12 0.01157108 0.00442433 2.615 0.009078 **
## my$bcidd1E2 0.00500682 0.00373187 1.342 0.180085
## my$bcidd1E3 0.00762603 0.00554298 1.376 0.169259
## my$bcidd1E4 0.00591987 0.00384425 1.540 0.123964
## my$bcidd1E5 -0.00657763 0.00498803 -1.319 0.187644
## my$bcidd1E6 -0.00139332 0.00372470 -0.374 0.708444
## my$bcidd1E7 0.00852322 0.00396829 2.148 0.032020 *
## my$bcidd1E8 -0.00369015 0.00443221 -0.833 0.405327
## my$bcidd1E9 0.00639603 0.00436339 1.466 0.143076
## my$bcidd1F1 -0.00012498 0.00381935 -0.033 0.973904
## my$bcidd1F10 0.00393912 0.00423405 0.930 0.352467
## my$bcidd1F11 0.00670256 0.00469224 1.428 0.153547
## my$bcidd1F12 0.00778245 0.00385357 2.020 0.043756 *
## my$bcidd1F2 0.00582870 0.00435141 1.339 0.180780
## my$bcidd1F3 0.00651265 0.00596554 1.092 0.275281
## my$bcidd1F4 -0.02625366 0.00923796 -2.842 0.004596 **
## my$bcidd1F5 -0.04834088 0.00429859 -11.246 < 0.0000000000000002 ***
## my$bcidd1F6 -0.00012930 0.00381509 -0.034 0.972972
## my$bcidd1F7 -0.01483489 0.00492987 -3.009 0.002700 **
## my$bcidd1F8 -0.00220117 0.00605837 -0.363 0.716455
## my$bcidd1F9 0.00029767 0.01166664 0.026 0.979651
## my$bcidd1G1 0.01006985 0.00654635 1.538 0.124376
## my$bcidd1G10 0.00690433 0.00358418 1.926 0.054408 .
## my$bcidd1G11 0.00299309 0.00557059 0.537 0.591205
## my$bcidd1G12 0.00147423 0.00344561 0.428 0.668867
## my$bcidd1G2 0.00615615 0.00569582 1.081 0.280096
## my$bcidd1G3 0.00097634 0.00390830 0.250 0.802795
## my$bcidd1G4 -0.00335673 0.00478773 -0.701 0.483433
## my$bcidd1G5 -0.01510445 0.00537164 -2.812 0.005043 **
## my$bcidd1G6 -0.00650477 0.00473025 -1.375 0.169463
## my$bcidd1G7 -0.00911754 0.00430618 -2.117 0.034534 *
## my$bcidd1G8 -0.00467010 0.00489124 -0.955 0.339967
## my$bcidd1G9 -0.00345997 0.00470652 -0.735 0.462463
## my$bcidd1H11 0.00188567 0.00387909 0.486 0.627018
## my$bcidd1H2 -0.00688077 0.00535513 -1.285 0.199193
## my$bcidd1H3 -0.00074760 0.00640698 -0.117 0.907138
## my$bcidd1H4 -0.03016766 0.00925002 -3.261 0.001155 **
## my$bcidd1H5 -0.00289320 0.00438300 -0.660 0.509377
## my$bcidd1H6 0.00117090 0.00365259 0.321 0.748621
## my$bcidd1H7 -0.00459377 0.00852914 -0.539 0.590312
## my$bcidd1H8 -0.00449698 0.00614172 -0.732 0.464255
## my$bcidd1H9 -0.00798745 0.00584794 -1.366 0.172359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.127 on 819 degrees of freedom
## Multiple R-squared: 0.2653, Adjusted R-squared: 0.1836
## F-statistic: 3.25 on 91 and 819 DF, p-value: < 0.00000000000000022
hist(residuals(mymod), breaks = 50) # residuals look good, check the other diagnostic plots too...
hist(expand_residuals(mymod, my$r), breaks = 50) # this more accurately represents the distribution of residuals in the model output.
plot(mymod) # view the diagnostic plots
paste0("Number of BCs with fitness different than pop mean = ", sum(p.adjust(coef(summary(mymod))[,
4], method = "fdr") < 0.05)) # number of barcodes with fit different than pop weighted mean.
## [1] "Number of BCs with fitness different than pop mean = 3"
paste0("RsMSE = ", sqrt(mean(residuals(mymod)^2))) # report rmse (requires unscalled lm model)
## [1] "RsMSE = 0.0175768866199638"
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S4.html")
| Â | my$rw corrected | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| d 1 A 10 | 0.00133 | -0.00595 – 0.00861 | 0.35687 | 0.72127977 |
| d 1 A 11 | -0.00762 | -0.01676 – 0.00151 | -1.63579 | 0.10226678 |
| d 1 A 2 | 0.00811 | 0.00035 – 0.01588 | 2.04902 | 0.04077841 |
| d 1 A 3 | 0.00494 | -0.00367 – 0.01356 | 1.12442 | 0.26116464 |
| d 1 A 4 | 0.00796 | -0.00084 – 0.01675 | 1.77364 | 0.07649506 |
| d 1 A 5 | 0.00103 | -0.00916 – 0.01123 | 0.19867 | 0.84257323 |
| d 1 A 6 | -0.00086 | -0.00787 – 0.00616 | -0.23972 | 0.81061012 |
| d 1 A 7 | 0.00344 | -0.00659 – 0.01348 | 0.67247 | 0.50147576 |
| d 1 A 8 | 0.00584 | -0.00229 – 0.01396 | 1.40710 | 0.15977841 |
| d 1 A 9 | 0.00724 | -0.00010 – 0.01459 | 1.93297 | 0.05358384 |
| d 1 B 1 | -0.00527 | -0.01512 – 0.00457 | -1.05048 | 0.29380661 |
| d 1 B 10 | -0.00821 | -0.02028 – 0.00386 | -1.33370 | 0.18267430 |
| d 1 B 11 | -0.01331 | -0.02333 – -0.00328 | -2.60193 | 0.00943744 |
| d 1 B 12 | -0.00742 | -0.01826 – 0.00343 | -1.34063 | 0.18041306 |
| d 1 B 2 | 0.00290 | -0.00556 – 0.01135 | 0.67171 | 0.50195688 |
| d 1 B 3 | 0.00354 | -0.00497 – 0.01206 | 0.81575 | 0.41487762 |
| d 1 B 4 | 0.00951 | 0.00125 – 0.01777 | 2.25665 | 0.02429251 |
| d 1 B 5 | 0.00962 | 0.00135 – 0.01788 | 2.28075 | 0.02281947 |
| d 1 B 6 | 0.00135 | -0.00849 – 0.01118 | 0.26840 | 0.78846186 |
| d 1 B 7 | -0.00398 | -0.01252 – 0.00455 | -0.91533 | 0.36028749 |
| d 1 B 8 | 0.00387 | -0.00467 – 0.01240 | 0.88790 | 0.37485596 |
| d 1 B 9 | -0.00197 | -0.01212 – 0.00817 | -0.38153 | 0.70290904 |
| d 1 C 1 | -0.00130 | -0.01157 – 0.00897 | -0.24854 | 0.80377774 |
| d 1 C 10 | -0.00201 | -0.01402 – 0.01001 | -0.32770 | 0.74322104 |
| d 1 C 11 | 0.00333 | -0.00675 – 0.01340 | 0.64716 | 0.51770702 |
| d 1 C 12 | -0.00496 | -0.01526 – 0.00535 | -0.94240 | 0.34626553 |
| d 1 C 2 | 0.00310 | -0.00554 – 0.01174 | 0.70303 | 0.48223454 |
| d 1 C 3 | 0.00171 | -0.00896 – 0.01237 | 0.31401 | 0.75359137 |
| d 1 C 4 | 0.00064 | -0.00808 – 0.00937 | 0.14452 | 0.88512380 |
| d 1 C 5 | -0.00293 | -0.01216 – 0.00630 | -0.62300 | 0.53345671 |
| d 1 C 6 | -0.00108 | -0.00894 – 0.00678 | -0.27021 | 0.78706863 |
| d 1 C 7 | -0.00224 | -0.01906 – 0.01457 | -0.26145 | 0.79381434 |
| d 1 C 8 | -0.00815 | -0.01866 – 0.00235 | -1.52132 | 0.12856573 |
| d 1 C 9 | 0.01983 | 0.00809 – 0.03157 | 3.31181 | 0.00096754 |
| d 1 D 1 | 0.00371 | -0.00443 – 0.01186 | 0.89388 | 0.37165067 |
| d 1 D 10 | 0.01036 | 0.00241 – 0.01832 | 2.55490 | 0.01080165 |
| d 1 D 11 | -0.00005 | -0.00752 – 0.00742 | -0.01308 | 0.98957011 |
| d 1 D 12 | 0.00302 | -0.00456 – 0.01060 | 0.78029 | 0.43544853 |
| d 1 D 2 | -0.01098 | -0.05482 – 0.03285 | -0.49113 | 0.62346368 |
| d 1 D 3 | -0.00440 | -0.01542 – 0.00662 | -0.78273 | 0.43401108 |
| d 1 D 4 | -0.00694 | -0.02138 – 0.00751 | -0.94144 | 0.34675951 |
| d 1 D 5 | -0.00327 | -0.01432 – 0.00778 | -0.57998 | 0.56208732 |
| d 1 D 6 | -0.00844 | -0.01930 – 0.00242 | -1.52337 | 0.12805285 |
| d 1 D 7 | -0.00168 | -0.01075 – 0.00738 | -0.36387 | 0.71605183 |
| d 1 D 8 | 0.00181 | -0.00598 – 0.00960 | 0.45468 | 0.64946246 |
| d 1 D 9 | 0.00014 | -0.01055 – 0.01084 | 0.02581 | 0.97941897 |
| d 1 E 1 | -0.01271 | -0.02297 – -0.00244 | -2.42696 | 0.01544074 |
| d 1 E 10 | -0.00148 | -0.00997 – 0.00702 | -0.34034 | 0.73369072 |
| d 1 E 11 | -0.00708 | -0.01504 – 0.00088 | -1.74357 | 0.08160859 |
| d 1 E 12 | 0.01157 | 0.00290 – 0.02024 | 2.61533 | 0.00907808 |
| d 1 E 2 | 0.00501 | -0.00231 – 0.01232 | 1.34164 | 0.18008470 |
| d 1 E 3 | 0.00763 | -0.00324 – 0.01849 | 1.37580 | 0.16925925 |
| d 1 E 4 | 0.00592 | -0.00161 – 0.01345 | 1.53993 | 0.12396449 |
| d 1 E 5 | -0.00658 | -0.01635 – 0.00320 | -1.31868 | 0.18764413 |
| d 1 E 6 | -0.00139 | -0.00869 – 0.00591 | -0.37408 | 0.70844388 |
| d 1 E 7 | 0.00852 | 0.00075 – 0.01630 | 2.14783 | 0.03201970 |
| d 1 E 8 | -0.00369 | -0.01238 – 0.00500 | -0.83258 | 0.40532715 |
| d 1 E 9 | 0.00640 | -0.00216 – 0.01495 | 1.46584 | 0.14307587 |
| d 1 F 1 | -0.00012 | -0.00761 – 0.00736 | -0.03272 | 0.97390392 |
| d 1 F 10 | 0.00394 | -0.00436 – 0.01224 | 0.93034 | 0.35246715 |
| d 1 F 11 | 0.00670 | -0.00249 – 0.01590 | 1.42844 | 0.15354746 |
| d 1 F 12 | 0.00778 | 0.00023 – 0.01534 | 2.01954 | 0.04375588 |
| d 1 F 2 | 0.00583 | -0.00270 – 0.01436 | 1.33950 | 0.18078013 |
| d 1 F 3 | 0.00651 | -0.00518 – 0.01820 | 1.09171 | 0.27528055 |
| d 1 F 4 | -0.02625 | -0.04436 – -0.00815 | -2.84193 | 0.00459550 |
| d 1 F 5 | -0.04834 | -0.05677 – -0.03992 | -11.24575 | <0.001 |
| d 1 F 6 | -0.00013 | -0.00761 – 0.00735 | -0.03389 | 0.97297191 |
| d 1 F 7 | -0.01483 | -0.02450 – -0.00517 | -3.00919 | 0.00269975 |
| d 1 F 8 | -0.00220 | -0.01408 – 0.00967 | -0.36333 | 0.71645463 |
| d 1 F 9 | 0.00030 | -0.02257 – 0.02316 | 0.02551 | 0.97965060 |
| d 1 G 1 | 0.01007 | -0.00276 – 0.02290 | 1.53824 | 0.12437600 |
| d 1 G 10 | 0.00690 | -0.00012 – 0.01393 | 1.92634 | 0.05440800 |
| d 1 G 11 | 0.00299 | -0.00793 – 0.01391 | 0.53730 | 0.59120517 |
| d 1 G 12 | 0.00147 | -0.00528 – 0.00823 | 0.42786 | 0.66886664 |
| d 1 G 2 | 0.00616 | -0.00501 – 0.01732 | 1.08082 | 0.28009589 |
| d 1 G 3 | 0.00098 | -0.00668 – 0.00864 | 0.24981 | 0.80279516 |
| d 1 G 4 | -0.00336 | -0.01274 – 0.00603 | -0.70111 | 0.48343252 |
| d 1 G 5 | -0.01510 | -0.02563 – -0.00458 | -2.81189 | 0.00504284 |
| d 1 G 6 | -0.00650 | -0.01578 – 0.00277 | -1.37514 | 0.16946284 |
| d 1 G 7 | -0.00912 | -0.01756 – -0.00068 | -2.11732 | 0.03453391 |
| d 1 G 8 | -0.00467 | -0.01426 – 0.00492 | -0.95479 | 0.33996652 |
| d 1 G 9 | -0.00346 | -0.01268 – 0.00576 | -0.73514 | 0.46246280 |
| d 1 H 11 | 0.00189 | -0.00572 – 0.00949 | 0.48611 | 0.62701817 |
| d 1 H 2 | -0.00688 | -0.01738 – 0.00362 | -1.28489 | 0.19919285 |
| d 1 H 3 | -0.00075 | -0.01331 – 0.01181 | -0.11669 | 0.90713761 |
| d 1 H 4 | -0.03017 | -0.04830 – -0.01204 | -3.26136 | 0.00115463 |
| d 1 H 5 | -0.00289 | -0.01148 – 0.00570 | -0.66010 | 0.50937740 |
| d 1 H 6 | 0.00117 | -0.00599 – 0.00833 | 0.32057 | 0.74862140 |
| d 1 H 7 | -0.00459 | -0.02131 – 0.01212 | -0.53860 | 0.59031160 |
| d 1 H 8 | -0.00450 | -0.01653 – 0.00754 | -0.73220 | 0.46425473 |
| d 1 H 9 | -0.00799 | -0.01945 – 0.00347 | -1.36586 | 0.17235858 |
| Observations | 910 | |||
| R2 / R2 adjusted | 0.265 / 0.184 | |||
# -----------------------------------------------------------------------------
# create temporary plotting dataframe
# -----------------------------------------
myres <- as.data.frame(coef(summary(mymod)), stringsAsFactors = F) # create a myres dataframe to aid in plotting.
myres$p.adjust <- p.adjust(myres$`Pr(>|t|)`, method = "fdr")
myres$bcid <- gsub("^.*?d", "", rownames(myres)) # fix barcode id names
myres$mrw <- ddply(my, ~bcid, function(my) weighted.mean(my$rw, my$r))[, 2]
myres$rmrw <- ddply(my, ~bcid, function(my) mean(my$r))[, 2]
myres$mcmrw <- myres$mrw - weighted.mean(myres$mrw, myres$rmrw)
myres$color <- "fdr > 0.05"
myres$color[myres$p.adjust <= 0.05] <- "fdr <= 0.05"
myres$color <- as.factor(myres$color)
myres$color <- relevel(myres$color, "fdr > 0.05")
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- myres$mcmrw
myr <- myres$rmrw
myc <- myres$color
mybinwidth <- 0.0025
mytitle <- NULL
myylab <- "Count"
myxlab <- "Fitness difference from population mean"
myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.3
mywidth <- 3.3
p <- ggplot(myres, aes(x = myx, color = myc, fill = myc)) + geom_histogram(alpha = myfillapha,
position = "identity", binwidth = mybinwidth) + geom_rug(alpha = myptalpha) +
xlim(c(-0.06, 0.06)) + ylim(c(NA, 15)) + geom_vline(xintercept = 0, linetype = "dashed",
color = myhighlight) + scale_color_manual(values = mypalette, guide = FALSE) +
scale_fill_manual(values = mypalette) + theme_classic() + theme(legend.position = c(0.8,
0.8)) + labs(fill = "Significance") + ggtitle(mytitle) + ylab(myylab) +
xlab(myxlab) + geom_hline(yintercept = -0.01, colour = "white", size = 1.01) +
theme(axis.title.y = element_text(size = mytextsize, face = "bold", margin = margin(0,
0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize,
angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) +
theme(legend.title = element_text(size = mytextsize, face = "bold")) + theme(legend.text = element_text(size = mytextsize))
p
## Warning: Removed 4 rows containing missing values (geom_bar).
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_2.pdf", width = mywidth, height = myheight)
p
## Warning: Removed 4 rows containing missing values (geom_bar).
dev.off() # one column wide, equal height
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(p, myres, mymod, my)
Analysis: Calculate sequenced library reads summary values
paste0("matching counts (match f and r primer barcode, and moby bc): ", sum(c(poc.ce,
poc.cc, faevo.1.ce, faevo.1.cc, faevo.2.ce, faevo.2.cc), na.rm = T)) # match forward primer, reverse primer, and yeast strain BCs
## [1] "matching counts (match f and r primer barcode, and moby bc): 104365740"
paste0("analysis counts: ", sum(c(poc.ce, poc.cc, faevo.ce, faevo.cc), na.rm = T)) # expected and cc counts in the analysis data
## [1] "analysis counts: 84776627"
paste0("expected analysis counts: ", sum(c(poc.ce, faevo.ce), na.rm = T)) # expected counts in the analysis data
## [1] "expected analysis counts: 81780343"
paste0("median counts per barcode: ", median(c(poc.ce, faevo.ce), na.rm = T)) # median number of counts per barcode in the expected counts dataset.
## [1] "median counts per barcode: 2961"
Analysis: Assess POC library barcode cross-contamination rate.
# calculations
# ----------------------------------------------------------------
cc = (((colSums(poc.cc, na.rm = T))/(colSums(poc.ce, na.rm = T) + colSums(poc.cc,
na.rm = T)))/colSums(!is.na(poc.cc))) * 100 # cc rate per exp
cc[is.nan(cc)] <- NA # NA nans
myt <- as.data.frame(cc) # convert to dataframe
myt$tc <- (colSums(poc.ce, na.rm = T) + colSums(poc.cc, na.rm = T)) # total counts to use as weight in weighted mean
cc001 <- weighted.mean(myt$cc, myt$tc, na.rm = T) # save weighted mean cc (by expected reads) for grand cc rate.
r001 <- mean(myt$tc) # save mean for below calculation of grand cc rate.
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- myt$cc
myr <- myt$tc
myc <- setpalette[1]
mybinwidth <- 0.01
mytitle <- NULL
myylab <- "Count"
myxlab <- "Mean Barcode Contamination Rate (Percent Total Counts)"
myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
p <- ggplot(myt, aes(x = myx)) + geom_histogram(alpha = myfillapha, color = myc,
fill = myc, binwidth = mybinwidth) + geom_vline(xintercept = weighted.mean(myx,
myr, na.rm = T), linetype = "dashed", color = myhighlight) + geom_rug(color = myc,
alpha = myptalpha) + xlim(NA, 1) + theme_classic() + ggtitle(mytitle) +
ylab(myylab) + xlab(myxlab) + theme(axis.title.y = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize,
angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) +
geom_hline(yintercept = -0.01, colour = "white", size = 1.01)
p
## Warning: Removed 20 rows containing non-finite values (stat_bin).
## Warning: Removed 1 rows containing missing values (geom_bar).
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S1_A.pdf", width = mywidth, height = myheight)
p
## Warning: Removed 20 rows containing non-finite values (stat_bin).
## Warning: Removed 1 rows containing missing values (geom_bar).
dev.off() # 2columns wide, half as tall
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(cc, myt, p) # remove the temporary variable
Analysis: Assess 250-generation experiment consensus library barcode cross-contamination rate.
# Calculations
# ----------------------------------------------------------------
cc = ((colSums(faevo.cc, na.rm = T))/(colSums(faevo.ce, na.rm = T) + colSums(faevo.cc,
na.rm = T)))/colSums(!is.na(faevo.cc)) * 100 # cc rate per exp
cc[is.nan(cc)] <- NA # NA nans
myt <- as.data.frame(cc) # convert to dataframe
myt$tc <- (colSums(faevo.ce, na.rm = T) + colSums(faevo.cc, na.rm = T)) # total counts to use as weight in weighted mean
cc002 <- weighted.mean(myt$cc, myt$tc, na.rm = T) # save weighted mean cc (by expected reads) for grand cc rate.
r002 <- mean(myt$tc) # save mean for below calculation of grand cc rate.
# overall cc rate
paste0("mean cross contam as % total counts = ", weighted.mean(c(cc001, cc002),
c(r001, r002)), " %") # weighted mean contam weighted.mean(c(contam rate poc, contam rate faevo), c(total counts poc, total counts faevo))
## [1] "mean cross contam as % total counts = 0.0427158666120247 %"
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- myt$cc
myr <- myt$tc
myc <- setpalette[1]
mybinwidth <- 0.01
mytitle <- NULL
myylab <- "Count"
myxlab <- "Mean Barcode Contamination Rate (Percent Total Counts)"
myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
p <- ggplot(myt, aes(x = myx)) + geom_histogram(alpha = myfillapha, color = myc,
fill = myc, binwidth = mybinwidth) + geom_vline(xintercept = weighted.mean(myx,
myr, na.rm = T), linetype = "dashed", color = myhighlight) + geom_rug(color = myc,
alpha = myptalpha) + xlim(NA, 1) + theme_classic() + ggtitle(mytitle) +
ylab(myylab) + xlab(myxlab) + theme(axis.title.y = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize,
angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) +
geom_hline(yintercept = -0.01, colour = "white", size = 1.01)
p
## Warning: Removed 1 rows containing missing values (geom_bar).
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S1_B.pdf", width = mywidth, height = myheight)
p
## Warning: Removed 1 rows containing missing values (geom_bar).
dev.off() # 2columns wide, half as tall
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(cc, myt, p, cc001, cc002, r001, r002) # remove the temporary variable
Analysis: Assess change in barcode cross-contamination rate for samples from the 250-generation experiment that were resequenced.
# Calculations & models
# -------------------------------------------------------
faevo.1.ce <- faevo.1.ce[, (colnames(faevo.1.ce) %in% colnames(faevo.2.ce))] # subset the run 1 data to only include those columns in the rerun dataset..expected counts..
faevo.1.cc <- faevo.1.cc[, (colnames(faevo.1.cc) %in% colnames(faevo.2.cc))] # cross contamination counts
faevo.1.ce <- faevo.1.ce[, colnames(faevo.1.ce)]
faevo.1.cc <- faevo.1.cc[, colnames(faevo.1.ce)] # reorder the run 1 data by column name
faevo.2.ce <- faevo.2.ce[, colnames(faevo.1.ce)]
faevo.2.cc <- faevo.2.cc[, colnames(faevo.1.ce)] # reorder the run 2 data by run 1 column name (for paired t test below)
myt1 <- ((colSums(faevo.1.cc, na.rm = T))/(colSums(faevo.1.ce, na.rm = T) +
colSums(faevo.1.cc, na.rm = T)))/colSums(!is.na(faevo.1.cc)) # cross contam run 1
myt2 <- ((colSums(faevo.2.cc, na.rm = T))/(colSums(faevo.2.ce, na.rm = T) +
colSums(faevo.2.cc, na.rm = T)))/colSums(!is.na(faevo.2.cc)) # cross contam run 2
myt1r <- (colSums(faevo.1.ce, na.rm = T) + colSums(faevo.1.cc, na.rm = T)) # reads run 1
myt2r <- (colSums(faevo.2.ce, na.rm = T) + colSums(faevo.2.cc, na.rm = T)) # reads run 2
# unweighted t-test results match weighted t-test qualitatively. see
# weighted t.test imemediately below. var.test(myt1, myt2) # variance is
# not equal mymod <- t.test(myt1, myt2, alternative = 'greater', paired =
# TRUE, var.equal = FALSE); mymod # t.test, paired, 1 sided (greater),
# unequal variance
myt <- myt1 - myt2 # calculatation for weighted t.test below (contam run 1 - contam run 2)
mytr <- 2/((1/myt1r) + (1/myt2r)) # reads for calculation below (harmonic mean of reads for run 1 and run 2)
mymod2 <- wtd.t.test(x = myt, y = 0, weight = mytr, alternative = "greater",
mean1 = T) # t.test, paired, 1 sided (greater), weighted
# -----------------------------------------------------------------------------
# create temporary plotting dataframe
# -----------------------------------------
myt1 <- as.data.frame(cbind(myt1, myt1r))
colnames(myt1) <- c("cc", "r")
myt1$seq.run <- 1 # run 1 df
myt2 <- as.data.frame(cbind(myt2, myt2r))
colnames(myt2) <- c("cc", "r")
myt2$seq.run <- 2 # run 2 df
myt <- rbind(myt1, myt2)
myt$seq.run <- as.factor(myt$seq.run) # combine
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- myt$seq.run
myy <- myt$cc
myr <- myt$r
myx1 <- myt1$cc
myr1 <- myt1$r
myx2 <- myt2$cc
myr2 <- myt2$r
mytitle <- NULL
myylab <- "Mean Barcode Contamination Rate (% Total Counts)"
myxlab <- "Sequenced Libary Run Number"
myoutline <- "grey50"
myfill <- "grey90"
myhighlight <- "grey50"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette
mytextsize <- 8
myheight <- 7
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(myt, aes(x = myx, y = myy, weight = myr, color = myx)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_segment(x = 0.75, xend = 1.25, y = weighted.mean(myx1,
myr1), yend = weighted.mean(myx1, myr1), color = myhighlight) + geom_segment(x = 1.75,
xend = 2.25, y = weighted.mean(myx2, myr2), yend = weighted.mean(myx2, myr2),
color = myhighlight) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha,
size = myptsize) + geom_rug(sides = "l", alpha = myptalpha) + scale_color_manual(values = mypalette) +
geom_hline(yintercept = 0, lty = "dotted") + theme_classic() + theme(legend.position = "none") +
ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + scale_x_discrete(breaks = c("1",
"2"), labels = c("Library 1", "Library 2")) + theme(axis.title.y = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize,
angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0))
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S2.pdf", width = mywidth, height = myheight)
p
dev.off() # 7by7
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(faevo.1.ce, faevo.1.cc, faevo.2.ce, faevo.2.cc, myt1, myt1r, myt2, myt2r,
myt, mytr, p) # clean up this section.
rm(faevo.cc, faevo.ce, poc.cc, poc.ce) # clean up for next major set of analyses (remove the counts datasets to keep our global environment clean)
Analysis: Assess whether detected fitness change for barcodes in the 250-generation experiment is affected by barcode-cross contamination rate
my <- myfa[order(myfa$bctID), ]
# model & diagnostic plots
# ----------------------------------------------------
mymod <- lmer(scale(deltafit) ~ scale(ccdeltafit) + (1 | treatment/bctID), weights = rdeltafit,
data = my, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
summary(mymod) # problems with convergence here...solved with an optimizer and some extra iterations.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(deltafit) ~ scale(ccdeltafit) + (1 | treatment/bctID)
## Data: my
## Weights: rdeltafit
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
##
## REML criterion at convergence: 797.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1167 -0.4622 0.0163 0.6366 3.5030
##
## Random effects:
## Groups Name Variance Std.Dev.
## bctID:treatment (Intercept) 0.1174 0.3426
## treatment (Intercept) 0.1865 0.4319
## Residual 368.2615 19.1901
## Number of obs: 608, groups: bctID:treatment, 152; treatment, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.11495 0.17940 4.94812 -0.641
## scale(ccdeltafit) -0.26885 0.01858 579.59718 -14.473
## Pr(>|t|)
## (Intercept) 0.55
## scale(ccdeltafit) <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scl(ccdltf) 0.007
# mymod <- lm(deltafit ~ ccdeltafit, weights = rdeltafit, data =
# my);r.squaredGLMM(mymod) # get the r2 for the lm version of the model if
# desired.
hist(residuals(mymod), breaks = 100) # residuals have outlying points -- but biew the weighted version of the residual histogram as well....
hist(expand_residuals(mymod, my$rdeltafit), breaks = 100) # this more accurately represents the distribution of residuals in the model output.
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
show.r2 = F, file = "000_Table_S6.html")
| Â | scale(deltafit) | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | -0.11495 | -0.46656 – 0.23667 | -0.64073 | 0.55018888 |
| scale(ccdeltafit) | -0.26885 | -0.30526 – -0.23244 | -14.47319 | <0.001 |
| Random Effects | ||||
| σ2 | 368.26 | |||
| τ00 bctID:treatment | 0.12 | |||
| τ00 treatment | 0.19 | |||
| ICC | 0.00 | |||
| N bctID | 152 | |||
| N treatment | 6 | |||
| Observations | 608 | |||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
mymod <- lmer(deltafit ~ ccdeltafit + (1 | treatment/bctID), weights = rdeltafit,
data = my, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))) # unmanipulated data for plotting
## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling
myx <- my$ccdeltafit
myy <- my$deltafit
myr <- my$rdeltafit
myc <- as.character(my$treatment)
myc[myc == "CM+Ethanoldiploid0.001"] <- "CM+EtOH \n Diploid \n 0.001 \n"
myc[myc == "CM+Saltdiploid0.001"] <- "CM+NaCl \n Diploid \n 0.001 \n"
myc[myc == "CMdiploid0.001"] <- "CM \n Diploid \n 0.001 \n"
myc[myc == "CMdiploid0.004"] <- "CM \n Diploid \n 0.004 \n"
myc[myc == "CMdiploid0.00025"] <- "CM \n Diploid \n 0.00025 \n"
myc[myc == "CMhaploid0.001"] <- "CM \n Haploid \n 0.001 \n"
myc <- factor(myc, levels = c("CM \n Diploid \n 0.001 \n", "CM+EtOH \n Diploid \n 0.001 \n",
"CM+NaCl \n Diploid \n 0.001 \n", "CM \n Diploid \n 0.00025 \n", "CM \n Diploid \n 0.004 \n",
"CM \n Haploid \n 0.001 \n"))
myg <- paste0(my$treatment, my$bctID)
mytitle <- NULL
myylab <- "Change in Fitness"
myxlab <- "Mean Barcode Contamination Rate (% Total Counts)"
myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette
mytextsize <- 8
myheight <- 7
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, color = myc)) + geom_line(aes(y = predict(mymod),
group = myg, color = myc), size = 1, alpha = myptalpha) + geom_point(aes(color = myc),
alpha = myptalpha, size = myptsize) + geom_rug(alpha = myptalpha) + scale_color_manual(values = setpalette) +
theme_classic() + labs(color = "Treatment") + ggtitle(mytitle) + ylab(myylab) +
xlab(myxlab) + theme(axis.title.y = element_text(size = mytextsize, face = "bold",
margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize,
angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) +
theme(legend.title = element_text(size = mytextsize, face = "bold")) + theme(legend.text = element_text(size = mytextsize))
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S3.pdf", width = mywidth, height = myheight)
p
dev.off() # 7by7
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p, weights)
Analysis: Assess fitness change in 250-generations of experimental evolution for barcodes in the 250-generation experiment.
my <- myfa[order(myfa$bctID), ]
# model & diagnostic plots
# ----------------------------------------------------
mymod <- lm(deltafit ~ bctID + 0, weights = rdeltafit, data = my)
summary(mymod)
##
## Call:
## lm(formula = deltafit ~ bctID + 0, data = my, weights = rdeltafit)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -8.9024 -0.9879 0.0734 0.9561 4.8229
##
## Coefficients:
## Estimate Std. Error t value
## bctIDCM+Ethanoldiploid0.001d1A2 0.00224666 0.01223756 0.184
## bctIDCM+Ethanoldiploid0.001d1A4 0.01765500 0.01365937 1.293
## bctIDCM+Ethanoldiploid0.001d1A5 -0.00183138 0.01601185 -0.114
## bctIDCM+Ethanoldiploid0.001d1A6 0.00847721 0.01160185 0.731
## bctIDCM+Ethanoldiploid0.001d1A7 0.00837054 0.01577339 0.531
## bctIDCM+Ethanoldiploid0.001d1D10 0.03354979 0.01054805 3.181
## bctIDCM+Ethanoldiploid0.001d1D11 0.02980938 0.01093362 2.726
## bctIDCM+Ethanoldiploid0.001d1D12 0.02399057 0.01018367 2.356
## bctIDCM+Ethanoldiploid0.001d1D9 0.02343344 0.02027118 1.156
## bctIDCM+Ethanoldiploid0.001d1E10 0.04807661 0.01120179 4.292
## bctIDCM+Ethanoldiploid0.001d1E11 0.03162314 0.01125755 2.809
## bctIDCM+Ethanoldiploid0.001d1E12 -0.03049231 0.01211023 -2.518
## bctIDCM+Ethanoldiploid0.001d1E9 -0.00999714 0.01070592 -0.934
## bctIDCM+Ethanoldiploid0.001d1F1 0.00938513 0.00879011 1.068
## bctIDCM+Ethanoldiploid0.001d1F2 0.01938389 0.01121524 1.728
## bctIDCM+Ethanoldiploid0.001d1G1 0.01657615 0.01500470 1.105
## bctIDCM+Ethanoldiploid0.001d1G2 -0.01127658 0.01818702 -0.620
## bctIDCM+Ethanoldiploid0.001d1H2 -0.01386016 0.01837215 -0.754
## bctIDCM+Ethanoldiploid0.001d1H4 0.00296372 0.02376768 0.125
## bctIDCM+Ethanoldiploid0.001d1H5 0.00347374 0.01381229 0.251
## bctIDCM+Ethanoldiploid0.001d1H6 0.00333766 0.01279047 0.261
## bctIDCM+Ethanoldiploid0.001d1H7 -0.00551949 0.02558074 -0.216
## bctIDCM+Saltdiploid0.001d1A2 0.04403600 0.01205900 3.652
## bctIDCM+Saltdiploid0.001d1A3 0.01686781 0.01180626 1.429
## bctIDCM+Saltdiploid0.001d1A4 0.02820342 0.01240104 2.274
## bctIDCM+Saltdiploid0.001d1A5 0.39852910 0.19667177 2.026
## bctIDCM+Saltdiploid0.001d1A7 0.02578701 0.01322947 1.949
## bctIDCM+Saltdiploid0.001d1F3 0.14079891 0.02033980 6.922
## bctIDCM+Saltdiploid0.001d1F4 0.18806044 0.02296431 8.189
## bctIDCM+Saltdiploid0.001d1F5 0.19913858 0.15799739 1.260
## bctIDCM+Saltdiploid0.001d1F6 0.23451078 0.01480773 15.837
## bctIDCM+Saltdiploid0.001d1F7 0.15255139 0.02121079 7.192
## bctIDCM+Saltdiploid0.001d1F8 0.12888704 0.01760927 7.319
## bctIDCM+Saltdiploid0.001d1G3 0.10101471 0.03918185 2.578
## bctIDCM+Saltdiploid0.001d1G4 0.13402619 0.04517263 2.967
## bctIDCM+Saltdiploid0.001d1G5 0.09585476 0.01936410 4.950
## bctIDCM+Saltdiploid0.001d1G6 0.08971402 0.03285475 2.731
## bctIDCM+Saltdiploid0.001d1G7 0.08740043 0.02106254 4.150
## bctIDCM+Saltdiploid0.001d1G8 0.12572385 0.03290001 3.821
## bctIDCM+Saltdiploid0.001d1H2 0.00495816 0.01539561 0.322
## bctIDCM+Saltdiploid0.001d1H3 0.02115633 0.01617586 1.308
## bctIDCM+Saltdiploid0.001d1H4 0.02137396 0.01932366 1.106
## bctIDCM+Saltdiploid0.001d1H5 0.00354951 0.01200179 0.296
## bctIDCM+Saltdiploid0.001d1H7 0.09417670 0.02624046 3.589
## bctIDCMdiploid0.00025d1A2 0.01559982 0.01599335 0.975
## bctIDCMdiploid0.00025d1A3 -0.00878718 0.01495776 -0.587
## bctIDCMdiploid0.00025d1A4 0.00174608 0.01441827 0.121
## bctIDCMdiploid0.00025d1A5 0.02677688 0.01511949 1.771
## bctIDCMdiploid0.00025d1A6 0.01198957 0.01303287 0.920
## bctIDCMdiploid0.00025d1A7 -0.00440271 0.01621492 -0.272
## bctIDCMdiploid0.00025d1B10 -0.00809898 0.01653299 -0.490
## bctIDCMdiploid0.00025d1B11 0.01723177 0.01242047 1.387
## bctIDCMdiploid0.00025d1B12 0.02983743 0.01830367 1.630
## bctIDCMdiploid0.00025d1B9 0.02155039 0.01336058 1.613
## bctIDCMdiploid0.00025d1C10 -0.01674094 0.01859239 -0.900
## bctIDCMdiploid0.00025d1C11 0.03216629 0.01622343 1.983
## bctIDCMdiploid0.00025d1C12 0.01496869 0.01326625 1.128
## bctIDCMdiploid0.00025d1C9 -0.00518047 0.01188006 -0.436
## bctIDCMdiploid0.00025d1D2 0.00370897 0.05550949 0.067
## bctIDCMdiploid0.00025d1E2 0.01289359 0.01052145 1.225
## bctIDCMdiploid0.00025d1H2 0.00760594 0.01762321 0.432
## bctIDCMdiploid0.00025d1H3 -0.00458899 0.01745253 -0.263
## bctIDCMdiploid0.00025d1H4 0.01843903 0.02426240 0.760
## bctIDCMdiploid0.00025d1H5 0.02541993 0.01377634 1.845
## bctIDCMdiploid0.00025d1H6 0.00414724 0.01352400 0.307
## bctIDCMdiploid0.00025d1H7 0.00190764 0.02782539 0.069
## bctIDCMdiploid0.001d1A11 0.03612225 0.01815435 1.990
## bctIDCMdiploid0.001d1A2 0.00809289 0.01417709 0.571
## bctIDCMdiploid0.001d1A2B 0.00927271 0.01214216 0.764
## bctIDCMdiploid0.001d1A3 0.00032174 0.01485237 0.022
## bctIDCMdiploid0.001d1A3B 0.00398055 0.01385717 0.287
## bctIDCMdiploid0.001d1A4 0.01780845 0.01310805 1.359
## bctIDCMdiploid0.001d1A4B 0.00153330 0.01338588 0.115
## bctIDCMdiploid0.001d1A5 0.02937872 0.01569201 1.872
## bctIDCMdiploid0.001d1A5B -0.00512256 0.01465871 -0.349
## bctIDCMdiploid0.001d1A6 0.02327336 0.01302785 1.786
## bctIDCMdiploid0.001d1A6B 0.00955064 0.01339667 0.713
## bctIDCMdiploid0.001d1A7 0.04290593 0.01344789 3.191
## bctIDCMdiploid0.001d1A8 0.01512575 0.01169949 1.293
## bctIDCMdiploid0.001d1A9 0.02268337 0.01128163 2.011
## bctIDCMdiploid0.001d1B1 0.05232713 0.01916650 2.730
## bctIDCMdiploid0.001d1B2 0.04063255 0.01363624 2.980
## bctIDCMdiploid0.001d1C1 0.05496896 0.01813523 3.031
## bctIDCMdiploid0.001d1C2 0.03217244 0.01299148 2.476
## bctIDCMdiploid0.001d1D4 0.08813997 0.05292837 1.665
## bctIDCMdiploid0.001d1D5 -0.01879019 0.01310839 -1.433
## bctIDCMdiploid0.001d1D6 0.04266468 0.04243456 1.005
## bctIDCMdiploid0.001d1D7 0.01221091 0.01202676 1.015
## bctIDCMdiploid0.001d1D8 -0.01270792 0.01081313 -1.175
## bctIDCMdiploid0.001d1E4 0.00301728 0.00782870 0.385
## bctIDCMdiploid0.001d1E5 0.03661495 0.02836404 1.291
## bctIDCMdiploid0.001d1E6 -0.02094510 0.00919326 -2.278
## bctIDCMdiploid0.001d1E7 0.02727065 0.01036588 2.631
## bctIDCMdiploid0.001d1E8 -0.00435115 0.01093102 -0.398
## bctIDCMdiploid0.001d1H11 0.03103102 0.01491899 2.080
## bctIDCMdiploid0.001d1H2 0.01207171 0.01797336 0.672
## bctIDCMdiploid0.001d1H2B -0.01771525 0.01722331 -1.029
## bctIDCMdiploid0.001d1H3 -0.00572467 0.01932551 -0.296
## bctIDCMdiploid0.001d1H3B 0.00004411 0.01824231 0.002
## bctIDCMdiploid0.001d1H4 0.00887392 0.02490203 0.356
## bctIDCMdiploid0.001d1H4B -0.00268735 0.02541084 -0.106
## bctIDCMdiploid0.001d1H5 0.01359003 0.01528546 0.889
## bctIDCMdiploid0.001d1H5B -0.00274906 0.01359686 -0.202
## bctIDCMdiploid0.001d1H6 0.03458717 0.01367055 2.530
## bctIDCMdiploid0.001d1H6B 0.02557412 0.01329414 1.924
## bctIDCMdiploid0.001d1H7 0.00942620 0.02647801 0.356
## bctIDCMdiploid0.001d1H8 0.01372130 0.02208970 0.621
## bctIDCMdiploid0.001d1H9 0.01922550 0.02317816 0.829
## bctIDCMdiploid0.004d1A2 0.00073845 0.01697648 0.043
## bctIDCMdiploid0.004d1A3 0.01219203 0.01903325 0.641
## bctIDCMdiploid0.004d1A4 0.05056340 0.01731766 2.920
## bctIDCMdiploid0.004d1A5 0.05149485 0.01798080 2.864
## bctIDCMdiploid0.004d1A6 0.01931859 0.01722893 1.121
## bctIDCMdiploid0.004d1A7 0.02865780 0.02004130 1.430
## bctIDCMdiploid0.004d1B3 0.02462782 0.01445493 1.704
## bctIDCMdiploid0.004d1B4 0.03125331 0.01659619 1.883
## bctIDCMdiploid0.004d1B5 0.04298949 0.01367318 3.144
## bctIDCMdiploid0.004d1B6 0.05327518 0.01738231 3.065
## bctIDCMdiploid0.004d1B7 0.07749931 0.01526527 5.077
## bctIDCMdiploid0.004d1C3 0.04273794 0.01469507 2.908
## bctIDCMdiploid0.004d1C4 0.02312320 0.01388803 1.665
## bctIDCMdiploid0.004d1C5 0.04334705 0.01293091 3.352
## bctIDCMdiploid0.004d1C6 0.06362223 0.01397915 4.551
## bctIDCMdiploid0.004d1C7 0.03835768 0.02896918 1.324
## bctIDCMdiploid0.004d1H2 0.02225147 0.02308198 0.964
## bctIDCMdiploid0.004d1H3 0.02106572 0.02173722 0.969
## bctIDCMdiploid0.004d1H4 0.00548235 0.03044485 0.180
## bctIDCMdiploid0.004d1H5 0.06026645 0.01911119 3.153
## bctIDCMdiploid0.004d1H6 0.04022208 0.01766259 2.277
## bctIDCMdiploid0.004d1H7 0.01164662 0.03155608 0.369
## bctIDCMhaploid0.001d1A2 0.17748276 0.18501465 0.959
## bctIDCMhaploid0.001d1A3 0.03475922 0.01642162 2.117
## bctIDCMhaploid0.001d1A3B 0.02653209 0.01137336 2.333
## bctIDCMhaploid0.001d1A4 0.07888998 0.01157424 6.816
## bctIDCMhaploid0.001d1A4B 0.00800659 0.01199416 0.668
## bctIDCMhaploid0.001d1A5 0.07567154 0.01371413 5.518
## bctIDCMhaploid0.001d1A5B 0.05376713 0.01259551 4.269
## bctIDCMhaploid0.001d1A6 0.02913551 0.01052784 2.767
## bctIDCMhaploid0.001d1A6B 0.02748147 0.00908708 3.024
## bctIDCMhaploid0.001d1A7 0.28424691 0.22971593 1.237
## bctIDCMhaploid0.001d1A7B 0.37549211 0.27832375 1.349
## bctIDCMhaploid0.001d1H2 0.05734718 0.01492184 3.843
## bctIDCMhaploid0.001d1H3 0.00588798 0.01683299 0.350
## bctIDCMhaploid0.001d1H3B 0.02271945 0.01846922 1.230
## bctIDCMhaploid0.001d1H4 0.06491278 0.03830520 1.695
## bctIDCMhaploid0.001d1H4B 0.01454690 0.01990335 0.731
## bctIDCMhaploid0.001d1H5 0.24925362 0.21028584 1.185
## bctIDCMhaploid0.001d1H5B 0.07768822 0.01659155 4.682
## bctIDCMhaploid0.001d1H6 0.10010351 0.02954796 3.388
## bctIDCMhaploid0.001d1H6B 0.11256605 0.02312772 4.867
## bctIDCMhaploid0.001d1H7 0.08865934 0.02084663 4.253
## bctIDCMhaploid0.001d1H7B 0.07939256 0.01932983 4.107
## Pr(>|t|)
## bctIDCM+Ethanoldiploid0.001d1A2 0.854419
## bctIDCM+Ethanoldiploid0.001d1A4 0.196832
## bctIDCM+Ethanoldiploid0.001d1A5 0.908990
## bctIDCM+Ethanoldiploid0.001d1A6 0.465352
## bctIDCM+Ethanoldiploid0.001d1A7 0.595903
## bctIDCM+Ethanoldiploid0.001d1D10 0.001570 **
## bctIDCM+Ethanoldiploid0.001d1D11 0.006650 **
## bctIDCM+Ethanoldiploid0.001d1D12 0.018906 *
## bctIDCM+Ethanoldiploid0.001d1D9 0.248287
## bctIDCM+Ethanoldiploid0.001d1E10 0.00002165288566661 ***
## bctIDCM+Ethanoldiploid0.001d1E11 0.005182 **
## bctIDCM+Ethanoldiploid0.001d1E12 0.012147 *
## bctIDCM+Ethanoldiploid0.001d1E9 0.350903
## bctIDCM+Ethanoldiploid0.001d1F1 0.286224
## bctIDCM+Ethanoldiploid0.001d1F2 0.084602 .
## bctIDCM+Ethanoldiploid0.001d1G1 0.269859
## bctIDCM+Ethanoldiploid0.001d1G2 0.535545
## bctIDCM+Ethanoldiploid0.001d1H2 0.450992
## bctIDCM+Ethanoldiploid0.001d1H4 0.900820
## bctIDCM+Ethanoldiploid0.001d1H5 0.801543
## bctIDCM+Ethanoldiploid0.001d1H6 0.794250
## bctIDCM+Ethanoldiploid0.001d1H7 0.829265
## bctIDCM+Saltdiploid0.001d1A2 0.000291 ***
## bctIDCM+Saltdiploid0.001d1A3 0.153770
## bctIDCM+Saltdiploid0.001d1A4 0.023413 *
## bctIDCM+Saltdiploid0.001d1A5 0.043309 *
## bctIDCM+Saltdiploid0.001d1A7 0.051883 .
## bctIDCM+Saltdiploid0.001d1F3 0.00000000001514729 ***
## bctIDCM+Saltdiploid0.001d1F4 0.00000000000000266 ***
## bctIDCM+Saltdiploid0.001d1F5 0.208173
## bctIDCM+Saltdiploid0.001d1F6 < 0.0000000000000002 ***
## bctIDCM+Saltdiploid0.001d1F7 0.00000000000263293 ***
## bctIDCM+Saltdiploid0.001d1F8 0.00000000000113446 ***
## bctIDCM+Saltdiploid0.001d1G3 0.010247 *
## bctIDCM+Saltdiploid0.001d1G4 0.003166 **
## bctIDCM+Saltdiploid0.001d1G5 0.00000104573824082 ***
## bctIDCM+Saltdiploid0.001d1G6 0.006566 **
## bctIDCM+Saltdiploid0.001d1G7 0.00003974952106726 ***
## bctIDCM+Saltdiploid0.001d1G8 0.000151 ***
## bctIDCM+Saltdiploid0.001d1H2 0.747562
## bctIDCM+Saltdiploid0.001d1H3 0.191568
## bctIDCM+Saltdiploid0.001d1H4 0.269265
## bctIDCM+Saltdiploid0.001d1H5 0.767557
## bctIDCM+Saltdiploid0.001d1H7 0.000368 ***
## bctIDCMdiploid0.00025d1A2 0.329882
## bctIDCMdiploid0.00025d1A3 0.557181
## bctIDCMdiploid0.00025d1A4 0.903664
## bctIDCMdiploid0.00025d1A5 0.077226 .
## bctIDCMdiploid0.00025d1A6 0.358086
## bctIDCMdiploid0.00025d1A7 0.786112
## bctIDCMdiploid0.00025d1B10 0.624463
## bctIDCMdiploid0.00025d1B11 0.166007
## bctIDCMdiploid0.00025d1B12 0.103764
## bctIDCMdiploid0.00025d1B9 0.107440
## bctIDCMdiploid0.00025d1C10 0.368373
## bctIDCMdiploid0.00025d1C11 0.048000 *
## bctIDCMdiploid0.00025d1C12 0.259774
## bctIDCMdiploid0.00025d1C9 0.662997
## bctIDCMdiploid0.00025d1D2 0.946757
## bctIDCMdiploid0.00025d1E2 0.221035
## bctIDCMdiploid0.00025d1H2 0.666246
## bctIDCMdiploid0.00025d1H3 0.792714
## bctIDCMdiploid0.00025d1H4 0.447657
## bctIDCMdiploid0.00025d1H5 0.065659 .
## bctIDCMdiploid0.00025d1H6 0.759244
## bctIDCMdiploid0.00025d1H7 0.945372
## bctIDCMdiploid0.001d1A11 0.047218 *
## bctIDCMdiploid0.001d1A2 0.568387
## bctIDCMdiploid0.001d1A2B 0.445453
## bctIDCMdiploid0.001d1A3 0.982727
## bctIDCMdiploid0.001d1A3B 0.774047
## bctIDCMdiploid0.001d1A4 0.174949
## bctIDCMdiploid0.001d1A4B 0.908856
## bctIDCMdiploid0.001d1A5 0.061817 .
## bctIDCMdiploid0.001d1A5B 0.726909
## bctIDCMdiploid0.001d1A6 0.074694 .
## bctIDCMdiploid0.001d1A6B 0.476265
## bctIDCMdiploid0.001d1A7 0.001518 **
## bctIDCMdiploid0.001d1A8 0.196716
## bctIDCMdiploid0.001d1A9 0.044951 *
## bctIDCMdiploid0.001d1B1 0.006576 **
## bctIDCMdiploid0.001d1B2 0.003039 **
## bctIDCMdiploid0.001d1C1 0.002576 **
## bctIDCMdiploid0.001d1C2 0.013632 *
## bctIDCMdiploid0.001d1D4 0.096546 .
## bctIDCMdiploid0.001d1D5 0.152415
## bctIDCMdiploid0.001d1D6 0.315226
## bctIDCMdiploid0.001d1D7 0.310496
## bctIDCMdiploid0.001d1D8 0.240516
## bctIDCMdiploid0.001d1E4 0.700111
## bctIDCMdiploid0.001d1E5 0.197395
## bctIDCMdiploid0.001d1E6 0.023170 *
## bctIDCMdiploid0.001d1E7 0.008807 **
## bctIDCMdiploid0.001d1E8 0.690775
## bctIDCMdiploid0.001d1H11 0.038086 *
## bctIDCMdiploid0.001d1H2 0.502150
## bctIDCMdiploid0.001d1H2B 0.304231
## bctIDCMdiploid0.001d1H3 0.767194
## bctIDCMdiploid0.001d1H3B 0.998072
## bctIDCMdiploid0.001d1H4 0.721741
## bctIDCMdiploid0.001d1H4B 0.915822
## bctIDCMdiploid0.001d1H5 0.374428
## bctIDCMdiploid0.001d1H5B 0.839864
## bctIDCMdiploid0.001d1H6 0.011740 *
## bctIDCMdiploid0.001d1H6B 0.055013 .
## bctIDCMdiploid0.001d1H7 0.722004
## bctIDCMdiploid0.001d1H8 0.534803
## bctIDCMdiploid0.001d1H9 0.407275
## bctIDCMdiploid0.004d1A2 0.965323
## bctIDCMdiploid0.004d1A3 0.522127
## bctIDCMdiploid0.004d1A4 0.003677 **
## bctIDCMdiploid0.004d1A5 0.004378 **
## bctIDCMdiploid0.004d1A6 0.262755
## bctIDCMdiploid0.004d1A7 0.153420
## bctIDCMdiploid0.004d1B3 0.089106 .
## bctIDCMdiploid0.004d1B4 0.060315 .
## bctIDCMdiploid0.004d1B5 0.001775 **
## bctIDCMdiploid0.004d1B6 0.002307 **
## bctIDCMdiploid0.004d1B7 0.00000056011850773 ***
## bctIDCMdiploid0.004d1C3 0.003811 **
## bctIDCMdiploid0.004d1C4 0.096605 .
## bctIDCMdiploid0.004d1C5 0.000868 ***
## bctIDCMdiploid0.004d1C6 0.00000684943418923 ***
## bctIDCMdiploid0.004d1C7 0.186138
## bctIDCMdiploid0.004d1H2 0.335547
## bctIDCMdiploid0.004d1H3 0.333005
## bctIDCMdiploid0.004d1H4 0.857174
## bctIDCMdiploid0.004d1H5 0.001720 **
## bctIDCMdiploid0.004d1H6 0.023234 *
## bctIDCMdiploid0.004d1H7 0.712242
## bctIDCMhaploid0.001d1A2 0.337921
## bctIDCMhaploid0.001d1A3 0.034829 *
## bctIDCMhaploid0.001d1A3B 0.020091 *
## bctIDCMhaploid0.001d1A4 0.00000000002976149 ***
## bctIDCMhaploid0.001d1A4B 0.504764
## bctIDCMhaploid0.001d1A5 0.00000005766073013 ***
## bctIDCMhaploid0.001d1A5B 0.00002392620167394 ***
## bctIDCMhaploid0.001d1A6 0.005879 **
## bctIDCMhaploid0.001d1A6B 0.002633 **
## bctIDCMhaploid0.001d1A7 0.216581
## bctIDCMhaploid0.001d1A7B 0.177968
## bctIDCMhaploid0.001d1H2 0.000139 ***
## bctIDCMhaploid0.001d1H3 0.726659
## bctIDCMhaploid0.001d1H3B 0.219285
## bctIDCMhaploid0.001d1H4 0.090830 .
## bctIDCMhaploid0.001d1H4B 0.465230
## bctIDCMhaploid0.001d1H5 0.236513
## bctIDCMhaploid0.001d1H5B 0.00000374585588773 ***
## bctIDCMhaploid0.001d1H6 0.000766 ***
## bctIDCMhaploid0.001d1H6B 0.00000156280464386 ***
## bctIDCMhaploid0.001d1H7 0.00002561158235811 ***
## bctIDCMhaploid0.001d1H7B 0.00004746151175424 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.922 on 456 degrees of freedom
## Multiple R-squared: 0.7186, Adjusted R-squared: 0.6248
## F-statistic: 7.66 on 152 and 456 DF, p-value: < 0.00000000000000022
hist(residuals(mymod), breaks = 100) # tails do look a little heavy
hist(expand_residuals(mymod, my$rdeltafit), breaks = 100) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # model fit (multi) -- looks pretty good.
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S7.html") # table does not have corrected p-values.
| Â | deltafit | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| CM+Ethanoldiploid 0.001 d 1 A 2 | 0.00225 | -0.02174 – 0.02623 | 0.18359 | 0.85441888 |
| CM+Ethanoldiploid 0.001 d 1 A 4 | 0.01765 | -0.00912 – 0.04443 | 1.29252 | 0.19683180 |
| CM+Ethanoldiploid 0.001 d 1 A 5 | -0.00183 | -0.03321 – 0.02955 | -0.11438 | 0.90898955 |
| CM+Ethanoldiploid 0.001 d 1 A 6 | 0.00848 | -0.01426 – 0.03122 | 0.73068 | 0.46535168 |
| CM+Ethanoldiploid 0.001 d 1 A 7 | 0.00837 | -0.02254 – 0.03929 | 0.53067 | 0.59590266 |
| CM+Ethanoldiploid 0.001 d 1 D 10 | 0.03355 | 0.01288 – 0.05422 | 3.18066 | 0.00156959 |
| CM+Ethanoldiploid 0.001 d 1 D 11 | 0.02981 | 0.00838 – 0.05124 | 2.72640 | 0.00664970 |
| CM+Ethanoldiploid 0.001 d 1 D 12 | 0.02399 | 0.00403 – 0.04395 | 2.35579 | 0.01890587 |
| CM+Ethanoldiploid 0.001 d 1 D 9 | 0.02343 | -0.01630 – 0.06316 | 1.15600 | 0.24828733 |
| CM+Ethanoldiploid 0.001 d 1 E 10 | 0.04808 | 0.02612 – 0.07003 | 4.29187 | 0.00002165 |
| CM+Ethanoldiploid 0.001 d 1 E 11 | 0.03162 | 0.00956 – 0.05369 | 2.80906 | 0.00518210 |
| CM+Ethanoldiploid 0.001 d 1 E 12 | -0.03049 | -0.05423 – -0.00676 | -2.51790 | 0.01214729 |
| CM+Ethanoldiploid 0.001 d 1 E 9 | -0.01000 | -0.03098 – 0.01099 | -0.93380 | 0.35090332 |
| CM+Ethanoldiploid 0.001 d 1 F 1 | 0.00939 | -0.00784 – 0.02661 | 1.06769 | 0.28622445 |
| CM+Ethanoldiploid 0.001 d 1 F 2 | 0.01938 | -0.00260 – 0.04137 | 1.72835 | 0.08460206 |
| CM+Ethanoldiploid 0.001 d 1 G 1 | 0.01658 | -0.01283 – 0.04598 | 1.10473 | 0.26985901 |
| CM+Ethanoldiploid 0.001 d 1 G 2 | -0.01128 | -0.04692 – 0.02437 | -0.62003 | 0.53554478 |
| CM+Ethanoldiploid 0.001 d 1 H 2 | -0.01386 | -0.04987 – 0.02215 | -0.75441 | 0.45099161 |
| CM+Ethanoldiploid 0.001 d 1 H 4 | 0.00296 | -0.04362 – 0.04955 | 0.12470 | 0.90081956 |
| CM+Ethanoldiploid 0.001 d 1 H 5 | 0.00347 | -0.02360 – 0.03055 | 0.25150 | 0.80154341 |
| CM+Ethanoldiploid 0.001 d 1 H 6 | 0.00334 | -0.02173 – 0.02841 | 0.26095 | 0.79424982 |
| CM+Ethanoldiploid 0.001 d 1 H 7 | -0.00552 | -0.05566 – 0.04462 | -0.21577 | 0.82926550 |
| CM+Saltdiploid 0.001 d 1 A 2 | 0.04404 | 0.02040 – 0.06767 | 3.65171 | 0.00029070 |
| CM+Saltdiploid 0.001 d 1 A 3 | 0.01687 | -0.00627 – 0.04001 | 1.42872 | 0.15376998 |
| CM+Saltdiploid 0.001 d 1 A 4 | 0.02820 | 0.00390 – 0.05251 | 2.27428 | 0.02341305 |
| CM+Saltdiploid 0.001 d 1 A 5 | 0.39853 | 0.01306 – 0.78400 | 2.02637 | 0.04330880 |
| CM+Saltdiploid 0.001 d 1 A 7 | 0.02579 | -0.00014 – 0.05172 | 1.94921 | 0.05188310 |
| CM+Saltdiploid 0.001 d 1 F 3 | 0.14080 | 0.10093 – 0.18066 | 6.92234 | <0.001 |
| CM+Saltdiploid 0.001 d 1 F 4 | 0.18806 | 0.14305 – 0.23307 | 8.18925 | <0.001 |
| CM+Saltdiploid 0.001 d 1 F 5 | 0.19914 | -0.11053 – 0.50881 | 1.26039 | 0.20817273 |
| CM+Saltdiploid 0.001 d 1 F 6 | 0.23451 | 0.20549 – 0.26353 | 15.83705 | <0.001 |
| CM+Saltdiploid 0.001 d 1 F 7 | 0.15255 | 0.11098 – 0.19412 | 7.19216 | <0.001 |
| CM+Saltdiploid 0.001 d 1 F 8 | 0.12889 | 0.09437 – 0.16340 | 7.31927 | <0.001 |
| CM+Saltdiploid 0.001 d 1 G 3 | 0.10101 | 0.02422 – 0.17781 | 2.57810 | 0.01024726 |
| CM+Saltdiploid 0.001 d 1 G 4 | 0.13403 | 0.04549 – 0.22256 | 2.96698 | 0.00316552 |
| CM+Saltdiploid 0.001 d 1 G 5 | 0.09585 | 0.05790 – 0.13381 | 4.95013 | 0.00000105 |
| CM+Saltdiploid 0.001 d 1 G 6 | 0.08971 | 0.02532 – 0.15411 | 2.73063 | 0.00656636 |
| CM+Saltdiploid 0.001 d 1 G 7 | 0.08740 | 0.04612 – 0.12868 | 4.14957 | 0.00003975 |
| CM+Saltdiploid 0.001 d 1 G 8 | 0.12572 | 0.06124 – 0.19021 | 3.82139 | 0.00015111 |
| CM+Saltdiploid 0.001 d 1 H 2 | 0.00496 | -0.02522 – 0.03513 | 0.32205 | 0.74756239 |
| CM+Saltdiploid 0.001 d 1 H 3 | 0.02116 | -0.01055 – 0.05286 | 1.30790 | 0.19156762 |
| CM+Saltdiploid 0.001 d 1 H 4 | 0.02137 | -0.01650 – 0.05925 | 1.10610 | 0.26926505 |
| CM+Saltdiploid 0.001 d 1 H 5 | 0.00355 | -0.01997 – 0.02707 | 0.29575 | 0.76755679 |
| CM+Saltdiploid 0.001 d 1 H 7 | 0.09418 | 0.04275 – 0.14561 | 3.58899 | 0.00036794 |
| C Mdiploid 0.00025 d 1 A 2 | 0.01560 | -0.01575 – 0.04695 | 0.97539 | 0.32988206 |
| C Mdiploid 0.00025 d 1 A 3 | -0.00879 | -0.03810 – 0.02053 | -0.58747 | 0.55718131 |
| C Mdiploid 0.00025 d 1 A 4 | 0.00175 | -0.02651 – 0.03001 | 0.12110 | 0.90366385 |
| C Mdiploid 0.00025 d 1 A 5 | 0.02678 | -0.00286 – 0.05641 | 1.77102 | 0.07722576 |
| C Mdiploid 0.00025 d 1 A 6 | 0.01199 | -0.01355 – 0.03753 | 0.91995 | 0.35808612 |
| C Mdiploid 0.00025 d 1 A 7 | -0.00440 | -0.03618 – 0.02738 | -0.27152 | 0.78611237 |
| C Mdiploid 0.00025 d 1 B 10 | -0.00810 | -0.04050 – 0.02431 | -0.48987 | 0.62446295 |
| C Mdiploid 0.00025 d 1 B 11 | 0.01723 | -0.00711 – 0.04158 | 1.38737 | 0.16600725 |
| C Mdiploid 0.00025 d 1 B 12 | 0.02984 | -0.00604 – 0.06571 | 1.63013 | 0.10376368 |
| C Mdiploid 0.00025 d 1 B 9 | 0.02155 | -0.00464 – 0.04774 | 1.61298 | 0.10744003 |
| C Mdiploid 0.00025 d 1 C 10 | -0.01674 | -0.05318 – 0.01970 | -0.90042 | 0.36837253 |
| C Mdiploid 0.00025 d 1 C 11 | 0.03217 | 0.00037 – 0.06396 | 1.98271 | 0.04799999 |
| C Mdiploid 0.00025 d 1 C 12 | 0.01497 | -0.01103 – 0.04097 | 1.12833 | 0.25977419 |
| C Mdiploid 0.00025 d 1 C 9 | -0.00518 | -0.02846 – 0.01810 | -0.43606 | 0.66299683 |
| C Mdiploid 0.00025 d 1 D 2 | 0.00371 | -0.10509 – 0.11251 | 0.06682 | 0.94675682 |
| C Mdiploid 0.00025 d 1 E 2 | 0.01289 | -0.00773 – 0.03352 | 1.22546 | 0.22103521 |
| C Mdiploid 0.00025 d 1 H 2 | 0.00761 | -0.02693 – 0.04215 | 0.43159 | 0.66624601 |
| C Mdiploid 0.00025 d 1 H 3 | -0.00459 | -0.03880 – 0.02962 | -0.26294 | 0.79271446 |
| C Mdiploid 0.00025 d 1 H 4 | 0.01844 | -0.02911 – 0.06599 | 0.75998 | 0.44765688 |
| C Mdiploid 0.00025 d 1 H 5 | 0.02542 | -0.00158 – 0.05242 | 1.84519 | 0.06565851 |
| C Mdiploid 0.00025 d 1 H 6 | 0.00415 | -0.02236 – 0.03065 | 0.30666 | 0.75924362 |
| C Mdiploid 0.00025 d 1 H 7 | 0.00191 | -0.05263 – 0.05644 | 0.06856 | 0.94537182 |
| C Mdiploid 0.001 d 1 A 11 | 0.03612 | 0.00054 – 0.07170 | 1.98973 | 0.04721759 |
| C Mdiploid 0.001 d 1 A 2 | 0.00809 | -0.01969 – 0.03588 | 0.57084 | 0.56838732 |
| C Mdiploid 0.001 d 1 A 2 B | 0.00927 | -0.01453 – 0.03307 | 0.76368 | 0.44545327 |
| C Mdiploid 0.001 d 1 A 3 | 0.00032 | -0.02879 – 0.02943 | 0.02166 | 0.98272673 |
| C Mdiploid 0.001 d 1 A 3 B | 0.00398 | -0.02318 – 0.03114 | 0.28726 | 0.77404701 |
| C Mdiploid 0.001 d 1 A 4 | 0.01781 | -0.00788 – 0.04350 | 1.35859 | 0.17494851 |
| C Mdiploid 0.001 d 1 A 4 B | 0.00153 | -0.02470 – 0.02777 | 0.11455 | 0.90885560 |
| C Mdiploid 0.001 d 1 A 5 | 0.02938 | -0.00138 – 0.06013 | 1.87221 | 0.06181742 |
| C Mdiploid 0.001 d 1 A 5 B | -0.00512 | -0.03385 – 0.02361 | -0.34945 | 0.72690929 |
| C Mdiploid 0.001 d 1 A 6 | 0.02327 | -0.00226 – 0.04881 | 1.78643 | 0.07469382 |
| C Mdiploid 0.001 d 1 A 6 B | 0.00955 | -0.01671 – 0.03581 | 0.71291 | 0.47626537 |
| C Mdiploid 0.001 d 1 A 7 | 0.04291 | 0.01655 – 0.06926 | 3.19053 | 0.00151809 |
| C Mdiploid 0.001 d 1 A 8 | 0.01513 | -0.00780 – 0.03806 | 1.29286 | 0.19671558 |
| C Mdiploid 0.001 d 1 A 9 | 0.02268 | 0.00057 – 0.04479 | 2.01065 | 0.04495102 |
| C Mdiploid 0.001 d 1 B 1 | 0.05233 | 0.01476 – 0.08989 | 2.73013 | 0.00657598 |
| C Mdiploid 0.001 d 1 B 2 | 0.04063 | 0.01391 – 0.06736 | 2.97975 | 0.00303894 |
| C Mdiploid 0.001 d 1 C 1 | 0.05497 | 0.01942 – 0.09051 | 3.03106 | 0.00257568 |
| C Mdiploid 0.001 d 1 C 2 | 0.03217 | 0.00671 – 0.05764 | 2.47643 | 0.01363235 |
| C Mdiploid 0.001 d 1 D 4 | 0.08814 | -0.01560 – 0.19188 | 1.66527 | 0.09654594 |
| C Mdiploid 0.001 d 1 D 5 | -0.01879 | -0.04448 – 0.00690 | -1.43345 | 0.15241540 |
| C Mdiploid 0.001 d 1 D 6 | 0.04266 | -0.04051 – 0.12583 | 1.00542 | 0.31522648 |
| C Mdiploid 0.001 d 1 D 7 | 0.01221 | -0.01136 – 0.03578 | 1.01531 | 0.31049575 |
| C Mdiploid 0.001 d 1 D 8 | -0.01271 | -0.03390 – 0.00849 | -1.17523 | 0.24051605 |
| C Mdiploid 0.001 d 1 E 4 | 0.00302 | -0.01233 – 0.01836 | 0.38541 | 0.70011148 |
| C Mdiploid 0.001 d 1 E 5 | 0.03661 | -0.01898 – 0.09221 | 1.29089 | 0.19739471 |
| C Mdiploid 0.001 d 1 E 6 | -0.02095 | -0.03896 – -0.00293 | -2.27831 | 0.02316980 |
| C Mdiploid 0.001 d 1 E 7 | 0.02727 | 0.00695 – 0.04759 | 2.63081 | 0.00880654 |
| C Mdiploid 0.001 d 1 E 8 | -0.00435 | -0.02578 – 0.01707 | -0.39806 | 0.69077548 |
| C Mdiploid 0.001 d 1 H 11 | 0.03103 | 0.00179 – 0.06027 | 2.07997 | 0.03808645 |
| C Mdiploid 0.001 d 1 H 2 | 0.01207 | -0.02316 – 0.04730 | 0.67164 | 0.50215015 |
| C Mdiploid 0.001 d 1 H 2 B | -0.01772 | -0.05147 – 0.01604 | -1.02856 | 0.30423060 |
| C Mdiploid 0.001 d 1 H 3 | -0.00572 | -0.04360 – 0.03215 | -0.29622 | 0.76719449 |
| C Mdiploid 0.001 d 1 H 3 B | 0.00004 | -0.03571 – 0.03580 | 0.00242 | 0.99807165 |
| C Mdiploid 0.001 d 1 H 4 | 0.00887 | -0.03993 – 0.05768 | 0.35635 | 0.72174082 |
| C Mdiploid 0.001 d 1 H 4 B | -0.00269 | -0.05249 – 0.04712 | -0.10576 | 0.91582228 |
| C Mdiploid 0.001 d 1 H 5 | 0.01359 | -0.01637 – 0.04355 | 0.88908 | 0.37442760 |
| C Mdiploid 0.001 d 1 H 5 B | -0.00275 | -0.02940 – 0.02390 | -0.20218 | 0.83986354 |
| C Mdiploid 0.001 d 1 H 6 | 0.03459 | 0.00779 – 0.06138 | 2.53005 | 0.01174027 |
| C Mdiploid 0.001 d 1 H 6 B | 0.02557 | -0.00048 – 0.05163 | 1.92371 | 0.05501270 |
| C Mdiploid 0.001 d 1 H 7 | 0.00943 | -0.04247 – 0.06132 | 0.35600 | 0.72200440 |
| C Mdiploid 0.001 d 1 H 8 | 0.01372 | -0.02957 – 0.05702 | 0.62116 | 0.53480274 |
| C Mdiploid 0.001 d 1 H 9 | 0.01923 | -0.02620 – 0.06465 | 0.82947 | 0.40727462 |
| C Mdiploid 0.004 d 1 A 2 | 0.00074 | -0.03253 – 0.03401 | 0.04350 | 0.96532350 |
| C Mdiploid 0.004 d 1 A 3 | 0.01219 | -0.02511 – 0.04950 | 0.64056 | 0.52212714 |
| C Mdiploid 0.004 d 1 A 4 | 0.05056 | 0.01662 – 0.08451 | 2.91976 | 0.00367654 |
| C Mdiploid 0.004 d 1 A 5 | 0.05149 | 0.01625 – 0.08674 | 2.86388 | 0.00437789 |
| C Mdiploid 0.004 d 1 A 6 | 0.01932 | -0.01445 – 0.05309 | 1.12129 | 0.26275545 |
| C Mdiploid 0.004 d 1 A 7 | 0.02866 | -0.01062 – 0.06794 | 1.42994 | 0.15341977 |
| C Mdiploid 0.004 d 1 B 3 | 0.02463 | -0.00370 – 0.05296 | 1.70377 | 0.08910621 |
| C Mdiploid 0.004 d 1 B 4 | 0.03125 | -0.00127 – 0.06378 | 1.88316 | 0.06031462 |
| C Mdiploid 0.004 d 1 B 5 | 0.04299 | 0.01619 – 0.06979 | 3.14407 | 0.00177490 |
| C Mdiploid 0.004 d 1 B 6 | 0.05328 | 0.01921 – 0.08734 | 3.06491 | 0.00230654 |
| C Mdiploid 0.004 d 1 B 7 | 0.07750 | 0.04758 – 0.10742 | 5.07684 | 0.00000056 |
| C Mdiploid 0.004 d 1 C 3 | 0.04274 | 0.01394 – 0.07154 | 2.90832 | 0.00381118 |
| C Mdiploid 0.004 d 1 C 4 | 0.02312 | -0.00410 – 0.05034 | 1.66497 | 0.09660494 |
| C Mdiploid 0.004 d 1 C 5 | 0.04335 | 0.01800 – 0.06869 | 3.35220 | 0.00086846 |
| C Mdiploid 0.004 d 1 C 6 | 0.06362 | 0.03622 – 0.09102 | 4.55122 | 0.00000685 |
| C Mdiploid 0.004 d 1 C 7 | 0.03836 | -0.01842 – 0.09514 | 1.32409 | 0.18613776 |
| C Mdiploid 0.004 d 1 H 2 | 0.02225 | -0.02299 – 0.06749 | 0.96402 | 0.33554719 |
| C Mdiploid 0.004 d 1 H 3 | 0.02107 | -0.02154 – 0.06367 | 0.96911 | 0.33300485 |
| C Mdiploid 0.004 d 1 H 4 | 0.00548 | -0.05419 – 0.06515 | 0.18007 | 0.85717391 |
| C Mdiploid 0.004 d 1 H 5 | 0.06027 | 0.02281 – 0.09772 | 3.15346 | 0.00171997 |
| C Mdiploid 0.004 d 1 H 6 | 0.04022 | 0.00560 – 0.07484 | 2.27725 | 0.02323376 |
| C Mdiploid 0.004 d 1 H 7 | 0.01165 | -0.05020 – 0.07350 | 0.36908 | 0.71224172 |
| C Mhaploid 0.001 d 1 A 2 | 0.17748 | -0.18514 – 0.54010 | 0.95929 | 0.33792087 |
| C Mhaploid 0.001 d 1 A 3 | 0.03476 | 0.00257 – 0.06694 | 2.11667 | 0.03482853 |
| C Mhaploid 0.001 d 1 A 3 B | 0.02653 | 0.00424 – 0.04882 | 2.33283 | 0.02009125 |
| C Mhaploid 0.001 d 1 A 4 | 0.07889 | 0.05620 – 0.10158 | 6.81600 | <0.001 |
| C Mhaploid 0.001 d 1 A 4 B | 0.00801 | -0.01550 – 0.03151 | 0.66754 | 0.50476447 |
| C Mhaploid 0.001 d 1 A 5 | 0.07567 | 0.04879 – 0.10255 | 5.51778 | 0.00000006 |
| C Mhaploid 0.001 d 1 A 5 B | 0.05377 | 0.02908 – 0.07845 | 4.26875 | 0.00002393 |
| C Mhaploid 0.001 d 1 A 6 | 0.02914 | 0.00850 – 0.04977 | 2.76747 | 0.00587909 |
| C Mhaploid 0.001 d 1 A 6 B | 0.02748 | 0.00967 – 0.04529 | 3.02424 | 0.00263330 |
| C Mhaploid 0.001 d 1 A 7 | 0.28425 | -0.16599 – 0.73448 | 1.23738 | 0.21658119 |
| C Mhaploid 0.001 d 1 A 7 B | 0.37549 | -0.17001 – 0.92100 | 1.34912 | 0.17796797 |
| C Mhaploid 0.001 d 1 H 2 | 0.05735 | 0.02810 – 0.08659 | 3.84317 | 0.00013869 |
| C Mhaploid 0.001 d 1 H 3 | 0.00589 | -0.02710 – 0.03888 | 0.34979 | 0.72665899 |
| C Mhaploid 0.001 d 1 H 3 B | 0.02272 | -0.01348 – 0.05892 | 1.23012 | 0.21928455 |
| C Mhaploid 0.001 d 1 H 4 | 0.06491 | -0.01016 – 0.13999 | 1.69462 | 0.09083022 |
| C Mhaploid 0.001 d 1 H 4 B | 0.01455 | -0.02446 – 0.05356 | 0.73088 | 0.46522961 |
| C Mhaploid 0.001 d 1 H 5 | 0.24925 | -0.16290 – 0.66141 | 1.18531 | 0.23651278 |
| C Mhaploid 0.001 d 1 H 5 B | 0.07769 | 0.04517 – 0.11021 | 4.68240 | 0.00000375 |
| C Mhaploid 0.001 d 1 H 6 | 0.10010 | 0.04219 – 0.15802 | 3.38783 | 0.00076555 |
| C Mhaploid 0.001 d 1 H 6 B | 0.11257 | 0.06724 – 0.15790 | 4.86715 | 0.00000156 |
| C Mhaploid 0.001 d 1 H 7 | 0.08866 | 0.04780 – 0.12952 | 4.25293 | 0.00002561 |
| C Mhaploid 0.001 d 1 H 7 B | 0.07939 | 0.04151 – 0.11728 | 4.10726 | 0.00004746 |
| Observations | 608 | |||
| R2 / R2 adjusted | 0.719 / 0.625 | |||
# -----------------------------------------------------------------------------
# create results dataframe
# ----------------------------------------------------
myres <- as.data.frame(coef(summary(mymod)), stringsAsFactors = F)
myres$p1s_l <- pt(coef(summary(mymod))[, 3], mymod$df, lower = T) # one-sided test p-values for lower
myres$p1s_h <- pt(coef(summary(mymod))[, 3], mymod$df, lower = F) # one-sided test p-values for higher
myres$p1s_l.adj <- p.adjust(myres$p1s_l, method = "fdr") # fdr adjusted p value (for 1-sided lower)
myres$p1s_h.adj <- p.adjust(myres$p1s_h, method = "fdr") # fdr adjusted p value (for 1-sided higher)
myres$bcid <- gsub("^.*?D", "", rownames(myres)) # fix barcode id names
myres$m_deltafit <- ddply(my, ~bctID, function(my) weighted.mean(my$deltafit,
my$rdeltafit))[, 2]
myres$m_rdeltafit <- ddply(my, ~bctID, function(my) mean(my$rdeltafit))[, 2]
myres$mcm_deltafit <- myres$m_deltafit - weighted.mean(myres$m_deltafit, myres$m_rdeltafit)
myres$color <- "fdr > 0.01"
myres$color[myres$p1s_h.adj <= 0.01] <- "fdr <= 0.01"
myres$color <- as.factor(myres$color)
myres$color <- relevel(myres$color, "fdr > 0.01")
myres$treatment <- gsub("0025d", "0025", gsub("001d", "001", gsub(".{4}$", "",
myres$bcid))) # strip bc id from treatment label, store in myres
# -----------------------------------------------------------------------------
# review results
# --------------------------------------------------------------
paste0("weighted mean fitness increase for increasers = ", weighted.mean(myres$m_deltafit[myres$p1s_h.adj <=
0.01], myres$m_rdeltafit[myres$p1s_h.adj <= 0.01])) # weighted mean fit inc of sig increasers only
## [1] "weighted mean fitness increase for increasers = 0.0679894955523455"
paste0("minimum fitness increase for increasers = ", min(myres$m_deltafit[myres$p1s_h.adj <=
0.01])) # minimum significant increase of 2.4%
## [1] "minimum fitness increase for increasers = 0.0274814710527647"
paste0("maximum fitness increase for increasers = ", max(myres$m_deltafit[myres$p1s_h.adj <=
0.01])) # maximum significant increase of 23.5%
## [1] "maximum fitness increase for increasers = 0.234510778858253"
# all treatments together
paste0("Number of barcodes with significant decreases in fitness = ", nrow(myres[myres$p1s_l.adj <=
0.01, ])) # check for significant fitness decrease: # find none
## [1] "Number of barcodes with significant decreases in fitness = 0"
paste0("Number of barcodes with significant increases in fitness = ", nrow(myres[myres$p1s_h.adj <=
0.01, ])) # check for significant fitness increase: # find 35 rows. (35 increasers out of 152)
## [1] "Number of barcodes with significant increases in fitness = 35"
paste0("Proportion of barcodes with significant increases in fitness = ", nrow(myres[myres$p1s_h.adj <=
0.01, ])/nrow(myres)) # prop increasers # 0.23
## [1] "Proportion of barcodes with significant increases in fitness = 0.230263157894737"
# CM diploid 1:1000
paste0("CM diploid 1:1000:")
## [1] "CM diploid 1:1000:"
myrest <- myres[myres$treatment == "CMdiploid0.001", ]
paste0("n barcodes in treatment = ", nrow(myrest))
## [1] "n barcodes in treatment = 42"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 3"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.0714285714285714"
# CM haploid 1:1000
paste0("CM haploid 1:1000:")
## [1] "CM haploid 1:1000:"
myrest <- myres[myres$treatment == "CMhaploid0.001", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 10"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.454545454545455"
# CM+salt diploid 1:1000
paste0("CM+salt diploid 1:1000:")
## [1] "CM+salt diploid 1:1000:"
myrest <- myres[myres$treatment == "CM+Saltdiploid0.001", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 11"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.5"
# CM+ethanol diploid 1:1000
paste0("CM+ethanol diploid 1:1000:")
## [1] "CM+ethanol diploid 1:1000:"
myrest <- myres[myres$treatment == "CM+Ethanoldiploid0.001", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 2"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.0909090909090909"
# CM diploid 1:250
paste0("CM diploid 1:250:")
## [1] "CM diploid 1:250:"
myrest <- myres[myres$treatment == "CMdiploid0.004", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 9"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.409090909090909"
# CM diploid 1:4000
paste0("CM diploid 1:4000:")
## [1] "CM diploid 1:4000:"
myrest <- myres[myres$treatment == "CMdiploid0.00025", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 0"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <=
0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0"
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- myres$m_deltafit
myr <- myres$m_rdeltafit
myc <- myres$color
mybinwidth <- 0.01
mytitle <- NULL
myylab <- "Count"
myxlab <- "Fitness Change"
myhighlight1 <- "black"
myhighlight2 <- setpalette[1]
myhighlight3 <- setpalette[2]
myfillapha <- 0.25
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.3
mywidth <- 3.3
p <- ggplot(myres, aes(x = myx, color = myc, fill = myc)) + geom_histogram(alpha = myfillapha,
position = "identity", binwidth = mybinwidth) + geom_vline(xintercept = weighted.mean(myx,
myr), linetype = "dashed", color = myhighlight) + geom_vline(xintercept = weighted.mean(myx[myc ==
"fdr > 0.01"], myr[myc == "fdr > 0.01"]), linetype = "dashed", color = myhighlight2) +
geom_vline(xintercept = weighted.mean(myx[myc == "fdr <= 0.01"], myr[myc ==
"fdr <= 0.01"]), linetype = "dashed", color = myhighlight3) + geom_rug(alpha = myptalpha) +
scale_color_manual(values = mypalette, guide = FALSE) + scale_fill_manual(values = mypalette) +
theme_classic() + theme(legend.position = c(0.8, 0.8)) + labs(fill = "Significance") +
ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + geom_hline(yintercept = -0.01,
colour = "white", size = 1.01) + theme(axis.title.y = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize,
angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) +
theme(legend.title = element_text(size = mytextsize, face = "bold")) + theme(legend.text = element_text(size = mytextsize))
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_3.pdf", width = mywidth, height = myheight)
p
dev.off() # 7by7
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, myres, myrest, p)
Analysis: Assess power to detect fitness change in this experimental evolution system, given replicate variance in the POC fitness assay dataset.
# create a power table for lms
# -------------------------------------------------
myrmse <- 1.7568866199638 # this is the rmse from the fitness equivalence model*100 --> effect size of 1 corresonds to an ~1.76 fitness change
myeffectsize <- 1/myrmse # effect size of 1 is a 1.76 fitnesss change, so for a fitness change of 1.00, expect effect size of 1/myrmse
mypower.lm <- as.data.frame(cbind(rep(NA, times = 12), NA, NA, NA, NA)) # create dataframe
colnames(mypower.lm) <- c("u", "v", "f2", "sig.level", "power") # fix column names
mypower.lm$u <- 1 # assign u
mypower.lm$v <- 6 # assign v
mypower.lm$sig.level <- c(0.05, 0.05, 0.05, 0.05, 0.01, 0.01, 0.01, 0.01, 0.005,
0.005, 0.005, 0.005) # assign significance
mypower.lm$f2 <- c(myeffectsize * 0.1, myeffectsize, myeffectsize * 2, myeffectsize *
5, myeffectsize * 0.1, myeffectsize, myeffectsize * 2, myeffectsize * 5,
myeffectsize * 0.1, myeffectsize, myeffectsize * 2, myeffectsize * 5) # Assign effect size
for (i in 1:nrow(mypower.lm)) {
# populate power column
mytemp <- pwr.f2.test(u = mypower.lm$u[i], v = mypower.lm$v[i], sig.level = mypower.lm$sig.level[i],
f2 = mypower.lm$f2[i])
mypower.lm$power[i] <- mytemp$power
}
# ------------------------------------------------------------------------------
# view and save
# ----------------------------------------------------------------
write.csv(mypower.lm, file = "000_Table_S5.csv")
# ------------------------------------------------------------------------------
rm(mytemp, i, mypower.lm)
Analysis: Assess effects of evoluitonary treatment on change in fitness in 250-generations of experimental evolution.
my <- myfa[order(myfa$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ----------------------------------------------------
mymod <- lmer(scale(deltafit) ~ scale(ccdeltafit) + treatment + (1 | bctID),
weights = rdeltafit, data = my, control = lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 20000)))
summary(mymod) # problems with convergence here...solved with an optimizer and some extra iterations.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(deltafit) ~ scale(ccdeltafit) + treatment + (1 | bctID)
## Data: my
## Weights: rdeltafit
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
##
## REML criterion at convergence: 789.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1056 -0.4647 0.0192 0.6301 3.5085
##
## Random effects:
## Groups Name Variance Std.Dev.
## bctID (Intercept) 0.1181 0.3437
## Residual 367.7128 19.1758
## Number of obs: 608, groups: bctID, 152
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.31415 0.05917 127.88336 -5.310
## scale(ccdeltafit) -0.27053 0.01860 576.86686 -14.548
## treatmentCM+Ethanoldiploid0.001 -0.29452 0.10043 127.81225 -2.932
## treatmentCM+Saltdiploid0.001 0.96163 0.10787 145.79857 8.915
## treatmentCMdiploid0.00025 -0.08752 0.10119 129.15660 -0.865
## treatmentCMdiploid0.004 0.34621 0.10249 138.20327 3.378
## treatmentCMhaploid0.001 0.27966 0.10888 134.62106 2.568
## Pr(>|t|)
## (Intercept) 0.00000047105690631 ***
## scale(ccdeltafit) < 0.0000000000000002 ***
## treatmentCM+Ethanoldiploid0.001 0.003987 **
## treatmentCM+Saltdiploid0.001 0.00000000000000184 ***
## treatmentCMdiploid0.00025 0.388701
## treatmentCMdiploid0.004 0.000949 ***
## treatmentCMhaploid0.001 0.011306 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(c) tCM+E0 tCM+S0 tCM0.000 tCM0.004
## scl(ccdltf) -0.020
## trCM+E0.001 -0.592 0.164
## trCM+S0.001 -0.546 -0.091 0.308
## trCM0.00025 -0.585 0.016 0.347 0.319
## trtmCM0.004 -0.576 -0.060 0.330 0.322 0.336
## trtmCM0.001 -0.545 0.112 0.338 0.288 0.319 0.307
# mymod <- lm(my$deltafit ~ my$ccdeltafit + my$treatment, weights =
# my$rdeltafit);r.squaredGLMM(mymod) # get the r2 for the lm version of the
# model if desired.
hist(residuals(mymod), breaks = 100)
hist(expand_residuals(mymod, my$rdeltafit), breaks = 100) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
anova(mymod) # test term significance
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------=
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
show.r2 = F, file = "000_Table_S8.html") # r2 on the lm version of htis model is 0.3059, adjusted = 0.299
| Â | scale(deltafit) | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | -0.31415 | -0.43012 – -0.19819 | -5.30957 | 0.00000047 |
| scale(ccdeltafit) | -0.27053 | -0.30697 – -0.23408 | -14.54794 | <0.001 |
| CM+Ethanoldiploid 0.001 | -0.29452 | -0.49136 – -0.09767 | -2.93240 | 0.00398690 |
| CM+Saltdiploid 0.001 | 0.96163 | 0.75021 – 1.17306 | 8.91470 | <0.001 |
| C Mdiploid 0.00025 | -0.08752 | -0.28584 – 0.11081 | -0.86489 | 0.38870062 |
| C Mdiploid 0.004 | 0.34621 | 0.14533 – 0.54710 | 3.37792 | 0.00094905 |
| C Mhaploid 0.001 | 0.27966 | 0.06625 – 0.49307 | 2.56846 | 0.01130558 |
| Random Effects | ||||
| σ2 | 367.71 | |||
| τ00 bctID | 0.12 | |||
| ICC | 0.00 | |||
| N bctID | 152 | |||
| Observations | 608 | |||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
mymod <- lmer((deltafit) ~ (ccdeltafit) + treatment + (1 | bctID), weights = rdeltafit,
data = my, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))) # problems with convergence here...solved with an optimizer and some extra iterations.
## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling
myx <- my$treatment
myy <- my$deltafit * 100
myr <- my$rdeltafit
myc <- my$treatment
mytitle <- NULL
myylab <- "Fitness Change"
myxlab <- "Treatment"
myoutline <- "grey20"
myfill <- "grey90"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() +
theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) +
scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001",
"CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004",
"CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text",
x = c(2, 3, 5, 6), y = c(14, 55, 15, 46), label = c("**", "***", "***",
"*"), size = 8, fontface = "bold")
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_4.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, weights)
Analysis: Calculate fixation rate summaries by treatment for the 250-generation experiment.
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# Review Results
# -------------------------------------------------------------- my[my$fixed
# == T,] # these are the fixed wells my$treatment[my$fixed == T] # these are
# the treatments of the fixed wells
paste0("proportion of wells that attained fixation (n=76 total wells)= ", nrow(my[my$fixed ==
T, ])/nrow(my)) # prop fixed
## [1] "proportion of wells that attained fixation (n=76 total wells)= 0.118421052631579"
# CM diploid 1:1000
myt <- my[my$treatment == "CMdiploid0.001", ]
paste0("n wells in CMdiploid0.001 = ", nrow(myt)) # wells
## [1] "n wells in CMdiploid0.001 = 21"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ])) # number fixed
## [1] "number of fixed wells in treatment = 1"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt)) # prop fixed
## [1] "proportion of fixed wells in treatment = 0.0476190476190476"
# CM haploid 1:1000
myt <- my[my$treatment == "CMhaploid0.001", ]
paste0("n wells in CMhaploid0.001 = ", nrow(myt)) # wells
## [1] "n wells in CMhaploid0.001 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ])) # number fixed
## [1] "number of fixed wells in treatment = 5"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt)) # prop fixed
## [1] "proportion of fixed wells in treatment = 0.454545454545455"
# CM+Salt diploid 1:1000
myt <- my[my$treatment == "CM+Saltdiploid0.001", ]
paste0("n wells in CM+Saltdiploid0.001 = ", nrow(myt)) # wells
## [1] "n wells in CM+Saltdiploid0.001 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ])) # number fixed
## [1] "number of fixed wells in treatment = 3"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt)) # prop fixed
## [1] "proportion of fixed wells in treatment = 0.272727272727273"
# CM+Ethanol 1:1000
myt <- my[my$treatment == "CM+Ethanoldiploid0.001", ]
paste0("n wells in CM+Ethanoldiploid0.001 = ", nrow(myt)) # wells
## [1] "n wells in CM+Ethanoldiploid0.001 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ])) # number fixed
## [1] "number of fixed wells in treatment = 0"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt)) # prop fixed
## [1] "proportion of fixed wells in treatment = 0"
# CM diploid 1:250
myt <- my[my$treatment == "CMdiploid0.004", ]
paste0("n wells in CMdiploid0.004 = ", nrow(myt)) # wells
## [1] "n wells in CMdiploid0.004 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ])) # number fixed
## [1] "number of fixed wells in treatment = 0"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt)) # prop fixed
## [1] "proportion of fixed wells in treatment = 0"
# CM diploid 1:4000
myt <- my[my$treatment == "CMdiploid0.00025", ]
paste0("n wells in CMdiploid0.00025 = ", nrow(myt)) # wells
## [1] "n wells in CMdiploid0.00025 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ])) # number fixed
## [1] "number of fixed wells in treatment = 0"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt)) # prop fixed
## [1] "proportion of fixed wells in treatment = 0"
# -----------------------------------------------------------------------------
rm(my, myt)
Analysis: Assess effects of evolutionary treatment on T-MAX (generation of maximum deviation from generation-0 barcoded strain abundance)
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tmaxdev ~ my$treatment
# + my$ccmaxdev_t0 + my$diff0, weights = my$rmaxdev_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$tmaxdev ~ my$treatment +
# my$diff0, weights = my$rmaxdev_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.
# without: mymod <- glm(my$tmaxdev ~ my$treatment + my$ccmaxdev, weights =
# my$rmaxdev); summary(mymod) # drop c as least significant term
mymod <- lm(my$tmaxdev ~ my$treatment, weights = my$rmaxdev)
summary(mymod) # FINAL MODEL!
##
## Call:
## lm(formula = my$tmaxdev ~ my$treatment, weights = my$rmaxdev)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -45513 -3616 1146 3910 16785
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 233.5992 14.7420 15.846
## my$treatmentCM+Ethanoldiploid0.001 -20.8306 17.1490 -1.215
## my$treatmentCM+Saltdiploid0.001 -48.7017 18.3948 -2.648
## my$treatmentCMdiploid0.00025 2.2063 20.8325 0.106
## my$treatmentCMdiploid0.004 -43.3530 27.8363 -1.557
## my$treatmentCMhaploid0.001 0.8156 16.9681 0.048
## Pr(>|t|)
## (Intercept) <0.0000000000000002 ***
## my$treatmentCM+Ethanoldiploid0.001 0.229
## my$treatmentCM+Saltdiploid0.001 0.010 *
## my$treatmentCMdiploid0.00025 0.916
## my$treatmentCMdiploid0.004 0.124
## my$treatmentCMhaploid0.001 0.962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9517 on 70 degrees of freedom
## Multiple R-squared: 0.1972, Adjusted R-squared: 0.1398
## F-statistic: 3.438 on 5 and 70 DF, p-value: 0.007779
anova(mymod) # test term significance with an lm
hist(residuals(mymod), breaks = 25)
hist(expand_residuals(mymod, my$rmaxdev), breaks = 25) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # check model fit
summary(mymod)$r.squared # for R2 value
## [1] 0.19716
summary(mymod)$adj.r.squared # for adjusted R2 value
## [1] 0.1398143
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S9.html")
| Â | my$tmaxdev | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | 233.59916 | 204.70539 – 262.49293 | 15.84584 | <0.001 |
| CM+Ethanoldiploid 0.001 | -20.83060 | -54.44209 – 12.78089 | -1.21468 | 0.22856965 |
| CM+Saltdiploid 0.001 | -48.70167 | -84.75481 – -12.64854 | -2.64758 | 0.01000874 |
| C Mdiploid 0.00025 | 2.20631 | -38.62458 – 43.03720 | 0.10591 | 0.91595897 |
| C Mdiploid 0.004 | -43.35303 | -97.91111 – 11.20506 | -1.55743 | 0.12387850 |
| C Mhaploid 0.001 | 0.81564 | -32.44114 – 34.07241 | 0.04807 | 0.96179819 |
| Observations | 76 | |||
| R2 / R2 adjusted | 0.197 / 0.140 | |||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tmaxdev
myr <- my$rmaxdev
myc <- my$treatment
mytitle <- NULL
myylab <- "T-MAX (Generation)"
myxlab <- "Treatment"
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
# geom_hline(yintercept = 0, lty = 'dotted') +
coord_flip() + theme_classic() + theme(legend.position = "none") + ggtitle(mytitle) +
ylab(myylab) + xlab(myxlab) + scale_x_discrete(labels = c("CM \n Diploid \n 0.001",
"CM+EtOH \n Diploid \n 0.001", "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025",
"CM \n Diploid \n 0.004", "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text",
x = c(3), y = c(255), label = "*", size = 8, fontface = "bold")
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S4.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)
Analysis: Assess effects of evolutionary treatment on M-MAX (magnitude of maximum deviation from generation-0 barcoded strain abundance)
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$mmaxdev ~ my$treatment
# + my$ccmaxdev_t0 + my$diff0, weights = my$rmaxdev_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$mmaxdev ~ my$treatment +
# my$diff0, weights = my$rmaxdev_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.
# without: mymod <- glm(my$mmaxdev ~ my$treatment + my$ccmaxdev, weights =
# my$rmaxdev); summary(mymod) # drop c as least significant term
mymod <- lm(my$mmaxdev ~ my$treatment, weights = my$rmaxdev)
summary(mymod) # FINAL MODEL!
##
## Call:
## lm(formula = my$mmaxdev ~ my$treatment, weights = my$rmaxdev)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -80.362 -11.841 -1.923 18.199 121.427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.22746 0.05437 4.184 0.0000820
## my$treatmentCM+Ethanoldiploid0.001 -0.06495 0.06325 -1.027 0.30799
## my$treatmentCM+Saltdiploid0.001 0.29240 0.06784 4.310 0.0000524
## my$treatmentCMdiploid0.00025 0.11117 0.07683 1.447 0.15240
## my$treatmentCMdiploid0.004 0.03904 0.10266 0.380 0.70491
## my$treatmentCMhaploid0.001 0.17791 0.06258 2.843 0.00585
##
## (Intercept) ***
## my$treatmentCM+Ethanoldiploid0.001
## my$treatmentCM+Saltdiploid0.001 ***
## my$treatmentCMdiploid0.00025
## my$treatmentCMdiploid0.004
## my$treatmentCMhaploid0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 35.1 on 70 degrees of freedom
## Multiple R-squared: 0.4564, Adjusted R-squared: 0.4176
## F-statistic: 11.76 on 5 and 70 DF, p-value: 0.00000002882
anova(mymod) # test term significance with an lm
hist(residuals(mymod), breaks = 25)
hist(expand_residuals(mymod, my$rmaxdev), breaks = 25) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # check model fit
summary(mymod)$r.squared # for R2 value
## [1] 0.4564297
summary(mymod)$adj.r.squared # for adjusted R2 value
## [1] 0.4176032
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S10.html")
| Â | my$mmaxdev | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | 0.22746 | 0.12089 – 0.33402 | 4.18350 | 0.00008197 |
| CM+Ethanoldiploid 0.001 | -0.06495 | -0.18891 – 0.05901 | -1.02693 | 0.30798766 |
| CM+Saltdiploid 0.001 | 0.29240 | 0.15943 – 0.42537 | 4.30998 | 0.00005236 |
| C Mdiploid 0.00025 | 0.11117 | -0.03942 – 0.26175 | 1.44685 | 0.15240115 |
| C Mdiploid 0.004 | 0.03904 | -0.16218 – 0.24025 | 0.38025 | 0.70491018 |
| C Mhaploid 0.001 | 0.17791 | 0.05526 – 0.30057 | 2.84299 | 0.00585276 |
| Observations | 76 | |||
| R2 / R2 adjusted | 0.456 / 0.418 | |||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$mmaxdev
myr <- my$rmaxdev
myc <- my$treatment
mytitle <- NULL
myylab <- "M-MAX"
myxlab <- "Treatment"
myoutline <- "grey20"
myfill <- "grey90"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() +
theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) +
scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001",
"CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004",
"CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text",
x = c(3, 6), y = c(0.91, 0.765), label = c("***", "**"), size = 8, fontface = "bold")
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S5.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)
Analysis: Assess effects of evolutionary treatment on T-MAX-RATE (generation of maximum rate of change in barcoded strain abundance)
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tmaxrate ~ my$treatment
# + my$ccmaxrate_t0 + my$diff0, weights = my$rmaxrate_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$tmaxrate ~ my$treatment +
# my$diff0, weights = my$rmaxrate_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.
# without: mymod <- glm(my$tmaxrate ~ my$treatment + my$ccmaxrate, weights =
# my$rmaxrate); summary(mymod) # drop c as least significant term
mymod <- lm(my$tmaxrate ~ my$treatment, weights = my$rmaxrate)
summary(mymod) # Final model!
##
## Call:
## lm(formula = my$tmaxrate ~ my$treatment, weights = my$rmaxrate)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -36500 -2877 985 3554 25478
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 235.251 7.709 30.518
## my$treatmentCM+Ethanoldiploid0.001 11.439 10.407 1.099
## my$treatmentCM+Saltdiploid0.001 -19.045 12.176 -1.564
## my$treatmentCMdiploid0.00025 9.910 17.741 0.559
## my$treatmentCMdiploid0.004 6.309 11.672 0.541
## my$treatmentCMhaploid0.001 -52.717 11.261 -4.681
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## my$treatmentCM+Ethanoldiploid0.001 0.275
## my$treatmentCM+Saltdiploid0.001 0.122
## my$treatmentCMdiploid0.00025 0.578
## my$treatmentCMdiploid0.004 0.591
## my$treatmentCMhaploid0.001 0.0000135 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9066 on 70 degrees of freedom
## Multiple R-squared: 0.3847, Adjusted R-squared: 0.3408
## F-statistic: 8.755 on 5 and 70 DF, p-value: 0.000001732
anova(mymod) # test term significance with an lm
hist(residuals(mymod), breaks = 25)
hist(expand_residuals(mymod, my$rmaxrate), breaks = 25) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # check model fit
summary(mymod)$r.squared # for R2 value
## [1] 0.3847431
summary(mymod)$adj.r.squared # for adjusted R2 value
## [1] 0.3407962
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S11.html")
| Â | my$tmaxrate | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | 235.25147 | 220.14297 – 250.35998 | 30.51820 | <0.001 |
| CM+Ethanoldiploid 0.001 | 11.43864 | -8.95850 – 31.83578 | 1.09914 | 0.27547276 |
| CM+Saltdiploid 0.001 | -19.04522 | -42.90918 – 4.81874 | -1.56420 | 0.12228104 |
| C Mdiploid 0.00025 | 9.90963 | -24.86202 – 44.68129 | 0.55857 | 0.57823532 |
| C Mdiploid 0.004 | 6.30932 | -16.56811 – 29.18674 | 0.54053 | 0.59054432 |
| C Mhaploid 0.001 | -52.71656 | -74.78791 – -30.64520 | -4.68130 | 0.00001354 |
| Observations | 76 | |||
| R2 / R2 adjusted | 0.385 / 0.341 | |||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tmaxrate
myr <- my$rmaxrate
myc <- my$treatment
mytitle <- NULL
myylab <- "T-MAX-RATE"
myxlab <- "Treatment"
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
# geom_hline(yintercept = 0, lty = 'dotted') +
coord_flip() + theme_classic() + theme(legend.position = "none") + ggtitle(mytitle) +
ylab(myylab) + xlab(myxlab) + scale_x_discrete(labels = c("CM \n Diploid \n 0.001",
"CM+EtOH \n Diploid \n 0.001", "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025",
"CM \n Diploid \n 0.004", "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text",
x = c(6), y = c(257), label = "***", size = 8, fontface = "bold")
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S6.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)
Analysis: Assess effects of evolutionary treatment on M-MAX-RATE (magnitude of maximum rate of change in barcoded strain abundance)
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$mmaxrate ~ my$treatment
# + my$ccmaxrate_t0 + my$diff0, weights = my$rmaxrate_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$mmaxrate ~ my$treatment +
# my$diff0, weights = my$rmaxrate_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.
# without: mymod <- glm(my$mmaxrate ~ my$treatment + my$ccmaxrate, weights =
# my$rmaxrate); summary(mymod) # drop c as least significant term
mymod <- lm(my$mmaxrate ~ my$treatment, weights = my$rmaxrate)
summary(mymod) # FINAL MODEL!
##
## Call:
## lm(formula = my$mmaxrate ~ my$treatment, weights = my$rmaxrate)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -4.0806 -0.8103 -0.1020 0.8827 4.2056
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.009582 0.001258 7.620
## my$treatmentCM+Ethanoldiploid0.001 -0.002035 0.001698 -1.199
## my$treatmentCM+Saltdiploid0.001 0.003154 0.001986 1.588
## my$treatmentCMdiploid0.00025 0.006305 0.002894 2.178
## my$treatmentCMdiploid0.004 -0.002101 0.001904 -1.104
## my$treatmentCMhaploid0.001 -0.002801 0.001837 -1.525
## Pr(>|t|)
## (Intercept) 0.0000000000916 ***
## my$treatmentCM+Ethanoldiploid0.001 0.2347
## my$treatmentCM+Saltdiploid0.001 0.1168
## my$treatmentCMdiploid0.00025 0.0328 *
## my$treatmentCMdiploid0.004 0.2736
## my$treatmentCMhaploid0.001 0.1319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.479 on 70 degrees of freedom
## Multiple R-squared: 0.2098, Adjusted R-squared: 0.1534
## F-statistic: 3.718 on 5 and 70 DF, p-value: 0.004837
anova(mymod) # test term significance with an lm
hist(residuals(mymod), breaks = 25)
hist(expand_residuals(mymod, my$rmaxrate), breaks = 25) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # check model fit
summary(mymod)$r.squared # for R2 value
## [1] 0.2098238
summary(mymod)$adj.r.squared # for adjusted R2 value
## [1] 0.1533827
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S12.html")
| Â | my$mmaxrate | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | 0.00958 | 0.00712 – 0.01205 | 7.61961 | <0.001 |
| CM+Ethanoldiploid 0.001 | -0.00204 | -0.00536 – 0.00129 | -1.19871 | 0.23468475 |
| CM+Saltdiploid 0.001 | 0.00315 | -0.00074 – 0.00705 | 1.58776 | 0.11684830 |
| C Mdiploid 0.00025 | 0.00630 | 0.00063 – 0.01198 | 2.17830 | 0.03275392 |
| C Mdiploid 0.004 | -0.00210 | -0.00583 – 0.00163 | -1.10353 | 0.27357441 |
| C Mhaploid 0.001 | -0.00280 | -0.00640 – 0.00080 | -1.52455 | 0.13187664 |
| Observations | 76 | |||
| R2 / R2 adjusted | 0.210 / 0.153 | |||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$mmaxrate
myr <- my$rmaxrate
myc <- my$treatment
mytitle <- NULL
myylab <- "M-MAX-RATE"
myxlab <- "Treatment"
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() +
theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) +
scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001",
"CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004",
"CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text",
x = c(4), y = c(0.0325), label = "*", size = 8, fontface = "bold")
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S7.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)
Analysis: Assess effects of evolutionary treatment on T-MAX-DIFF (generation of maximum difference in sympatric barcoded strain abundance)
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tmaxdiff ~ my$treatment
# + my$ccmaxdiff_t0 + my$diff0, weights = my$rmaxdiff_rt0); summary(mymod) #
# drop c as least significant term
mymod <- lm(my$tmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
summary(mymod) # # FINAL MODEL!
##
## Call:
## lm(formula = my$tmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -58426 -5302 574 8332 59655
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 241.970 34.101 7.096
## my$treatmentCM+Ethanoldiploid0.001 -22.732 35.785 -0.635
## my$treatmentCM+Saltdiploid0.001 -5.533 38.822 -0.143
## my$treatmentCMdiploid0.00025 11.312 77.413 0.146
## my$treatmentCMdiploid0.004 -23.716 63.676 -0.372
## my$treatmentCMhaploid0.001 1.612 35.479 0.045
## my$diff0 -328.617 81.224 -4.046
## Pr(>|t|)
## (Intercept) 0.000000000892 ***
## my$treatmentCM+Ethanoldiploid0.001 0.527381
## my$treatmentCM+Saltdiploid0.001 0.887074
## my$treatmentCMdiploid0.00025 0.884245
## my$treatmentCMdiploid0.004 0.710705
## my$treatmentCMhaploid0.001 0.963902
## my$diff0 0.000134 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19160 on 69 degrees of freedom
## Multiple R-squared: 0.2342, Adjusted R-squared: 0.1676
## F-statistic: 3.517 on 6 and 69 DF, p-value: 0.004283
anova(mymod) # test term significance with an lm
hist(residuals(mymod), breaks = 25)
hist(expand_residuals(mymod, my$rtcc), breaks = 25) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # check model fit
summary(mymod)$r.squared # for R2 value
## [1] 0.2342223
summary(mymod)$adj.r.squared # for adjusted R2 value
## [1] 0.1676329
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S13.html")
|  | my\(tmaxdiff</th> </tr> <tr> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; text-align:left; ">Predictors</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Estimates</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">CI</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Statistic</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">p</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">(Intercept)</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">241.97046</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">175.13346 – 308.80745</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">7.09567</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong><0.001</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Ethanoldiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-22.73154</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-92.86866 – 47.40558</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.63523</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.52738134</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Saltdiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-5.53349</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-81.62393 – 70.55695</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.14253</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.88707398</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.00025</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">11.31234</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-140.41438 – 163.03905</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.14613</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.88424511</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.004</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-23.71583</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-148.51929 – 101.08763</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.37244</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.71070458</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mhaploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">1.61153</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-67.92574 – 71.14880</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.04542</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.96390191</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">my\)diff 0 | -328.61667 | -487.81355 – -169.41979 | -4.04579 | 0.00013406 | |||
|---|---|---|---|---|---|---|---|---|
| Observations | 76 | |||||||
| R2 / R2 adjusted | 0.234 / 0.168 | |||||||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tmaxdiff
myr <- my$rmaxdiff_rt0
myc <- my$treatment
mytitle <- NULL
myylab <- "T-MAX-DIFF"
myxlab <- "Treatment"
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() +
theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) +
scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001",
"CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004",
"CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0))
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S8.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)
Analysis: Assess effects of evolutionary treatment on M-MAX-DIFF (magnitude of maximum difference in sympatric barcoded strain abundance)
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$mmaxdiff ~ my$treatment
# + my$ccmaxdiff_t0 + my$diff0, weights = my$rmaxdiff_rt0); summary(mymod) #
# drop c as least significant term
mymod <- lm(my$mmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
summary(mymod) # # FINAL MODEL!
##
## Call:
## lm(formula = my$mmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -90.008 -8.179 2.215 10.833 47.792
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.18985 0.04004 4.741 0.00001105
## my$treatmentCM+Ethanoldiploid0.001 0.00454 0.04202 0.108 0.914271
## my$treatmentCM+Saltdiploid0.001 0.17472 0.04559 3.833 0.000277
## my$treatmentCMdiploid0.00025 0.15524 0.09090 1.708 0.092165
## my$treatmentCMdiploid0.004 0.03311 0.07477 0.443 0.659311
## my$treatmentCMhaploid0.001 0.20429 0.04166 4.904 0.00000601
## my$diff0 0.29644 0.09538 3.108 0.002734
##
## (Intercept) ***
## my$treatmentCM+Ethanoldiploid0.001
## my$treatmentCM+Saltdiploid0.001 ***
## my$treatmentCMdiploid0.00025 .
## my$treatmentCMdiploid0.004
## my$treatmentCMhaploid0.001 ***
## my$diff0 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.49 on 69 degrees of freedom
## Multiple R-squared: 0.4789, Adjusted R-squared: 0.4336
## F-statistic: 10.57 on 6 and 69 DF, p-value: 0.00000002701
anova(mymod) # test term significance with an lm
hist(residuals(mymod), breaks = 25)
hist(expand_residuals(mymod, my$rtcc), breaks = 25) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # check model fit
summary(mymod)$r.squared # for R2 value
## [1] 0.4789295
summary(mymod)$adj.r.squared # for adjusted R2 value
## [1] 0.433619
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S14.html")
|  | my\(mmaxdiff</th> </tr> <tr> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; text-align:left; ">Predictors</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Estimates</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">CI</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Statistic</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">p</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">(Intercept)</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.18985</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.11137 – 0.26833</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">4.74128</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>0.00001105</strong></td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Ethanoldiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.00454</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.07782 – 0.08690</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.10805</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.91427138</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Saltdiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.17472</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.08537 – 0.26406</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">3.83269</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>0.00027655</strong></td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.00025</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.15524</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.02292 – 0.33340</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">1.70782</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.09216521</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.004</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.03311</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.11344 – 0.17965</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.44278</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.65931134</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mhaploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.20429</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.12264 – 0.28594</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">4.90375</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>0.00000601</strong></td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">my\)diff 0 | 0.29644 | 0.10951 – 0.48337 | 3.10815 | 0.00273414 | |||
|---|---|---|---|---|---|---|---|---|
| Observations | 76 | |||||||
| R2 / R2 adjusted | 0.479 / 0.434 | |||||||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$mmaxdiff
myr <- my$rmaxdiff_rt0
myc <- my$treatment
mytitle <- NULL
myylab <- "M-MAX-DIFF"
myxlab <- "Treatment"
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() +
theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) +
scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001",
"CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004",
"CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text",
x = c(3, 6), y = c(0.525), label = "***", size = 8, fontface = "bold")
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S9.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)
Analysis: Assess effects of evolutionary treatment on total change in barcode abundance (summed across six intervals)
my <- myevo[order(myevo$bctID), ] # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tcc ~ my$treatment +
# my$cctcc_t0 + my$diff0, weights = my$rtcc_rt0); summary(mymod) # drop c as
# least significant term mymod <- glm(my$tcc ~ my$treatment + my$diff0,
# weights = my$rtcc_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.
# without: mymod <- glm(my$tcc ~ my$treatment + my$cctcc, weights =
# my$rtcc); summary(mymod) # drop c as least significant term
mymod <- lm(my$tcc ~ my$treatment, weights = my$rtcc)
summary(mymod) # FINAL MODEL!
##
## Call:
## lm(formula = my$tcc ~ my$treatment, weights = my$rtcc)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -5.7198 -1.3007 -0.0058 1.0657 5.4905
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.0219194 0.0027409 7.997
## my$treatmentCM+Ethanoldiploid0.001 -0.0072401 0.0032633 -2.219
## my$treatmentCM+Saltdiploid0.001 0.0046752 0.0035042 1.334
## my$treatmentCMdiploid0.00025 0.0182730 0.0056216 3.251
## my$treatmentCMdiploid0.004 -0.0007325 0.0043115 -0.170
## my$treatmentCMhaploid0.001 -0.0066256 0.0033076 -2.003
## Pr(>|t|)
## (Intercept) 0.0000000000185 ***
## my$treatmentCM+Ethanoldiploid0.001 0.02976 *
## my$treatmentCM+Saltdiploid0.001 0.18647
## my$treatmentCMdiploid0.00025 0.00177 **
## my$treatmentCMdiploid0.004 0.86559
## my$treatmentCMhaploid0.001 0.04904 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.15 on 70 degrees of freedom
## Multiple R-squared: 0.3746, Adjusted R-squared: 0.33
## F-statistic: 8.387 on 5 and 70 DF, p-value: 0.000002951
anova(mymod) # test term significance with an lm
hist(residuals(mymod), breaks = 25)
hist(expand_residuals(mymod, my$rtcc), breaks = 25) # this more accurately represents the distribution of residuals in the model output. -- tails look fine here.
plot(mymod) # check model fit
summary(mymod)$r.squared # for R2 value
## [1] 0.3746437
summary(mymod)$adj.r.squared # for adjusted R2 value
## [1] 0.3299754
# -----------------------------------------------------------------------------
# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100,
file = "000_Table_S15.html")
| Â | my$tcc | |||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | 0.02192 | 0.01655 – 0.02729 | 7.99719 | <0.001 |
| CM+Ethanoldiploid 0.001 | -0.00724 | -0.01364 – -0.00084 | -2.21861 | 0.02975813 |
| CM+Saltdiploid 0.001 | 0.00468 | -0.00219 – 0.01154 | 1.33418 | 0.18646947 |
| C Mdiploid 0.00025 | 0.01827 | 0.00725 – 0.02929 | 3.25053 | 0.00177244 |
| C Mdiploid 0.004 | -0.00073 | -0.00918 – 0.00772 | -0.16989 | 0.86559011 |
| C Mhaploid 0.001 | -0.00663 | -0.01311 – -0.00014 | -2.00311 | 0.04903884 |
| Observations | 76 | |||
| R2 / R2 adjusted | 0.375 / 0.330 | |||
# -----------------------------------------------------------------------------
# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tcc
myr <- my$rtcc
myc <- my$treatment
mytitle <- NULL
myylab <- "Total Change in Abundance"
myxlab <- "Treatment"
mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) *
10) * ((myheight * mywidth)/(7^2)))) + 1
p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline,
fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy),
trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) +
geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"],
myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9,
xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]),
color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx ==
"CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) +
geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"],
myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) +
geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"],
myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx ==
"CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) +
geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"],
myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b",
position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) +
geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() +
theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) +
scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001",
"CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004",
"CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8,
face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8,
angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text",
x = c(2, 4, 6), y = c(0.0275, 0.06825, 0.026), label = c("*", "**", "*"),
size = 8, fontface = "bold")
p
# -----------------------------------------------------------------------------
# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S10.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen
## 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)